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Introduction 


An  important  initial  screening  step  in  the  detection  of  breast  caneer  is  the  ability  to  identify  select 
areas  of  atypieal  density  that  require  further  evaluation.  Currently,  mammography  is  the  elinieal 
standard  for  screening  and  provides  useful  but  at  times  ambiguous  information,  whieh  can 
necessitate  further  invasive  workup  of  benign  lesions.  Alternative  methods  sueh  as  elastography 
have  shown  potential  in  enhaneing  the  diagnostie  process  by  providing  information  about  the 
tissue  eomposition  [1,2].  Modality-independent  elastography  (MIE)  is  a  novel  image  processing 
technique  that  eombines  finite  element  models  of  soft-tissue  deformation  with  measures  of  image 
similarity  in  order  to  reconstruct  elastic  property  distributions  throughout  the  tissue.  The  basie 
requirements  for  the  method  are  two  images  of  the  tissue  in  different  states  of  deformation  (e.g. 
compression).  MIE  then  updates  the  estimate  of  the  material  properties  via  a  matehing  process 
between  the  two  images.  The  final  result  is  a  map  of  the  breast  (or  other  tissue  of  interest)  that 
refleets  material  inhomogeneity,  such  as  in  the  case  of  a  tumor  mass  that  disrupts  the  surrounding 
structure  of  normal  tissue.  Because  MIE  works  on  probing  the  differenees  between  images,  it  ean 
be  used  to  not  only  work  in  eoncert  with  more  traditional  screening  techniques  but  also  address  a 
possible  gap  when  those  methods  are  unable  to  directly  discern  tissues  of  interest. 


Body 

As  stated  in  the  original  proposal,  three  main  aims  of  this  projeet  are  to  (1)  expand  and  refine  the 
current  MIE  technique  to  enhanee  its  efficiency  and  capabilities,  (2)  to  perform  analyses  on 
texture  in  input  images  and  quantify  statistieal  parameters  capable  of  estimating  and  evaluating 
the  success  of  elastographie  reconstruction,  and  (3)  to  engineer  a  deviee  that  can  accurately 
produee  eompressive  forees  necessary  for  phantom  setups  within  eurrent  imaging  systems, 
providing  the  basis  for  a  future  deviee  that  can  be  used  in  a  elinieal  setting.  In  this  past  year, 
progress  on  all  three  aims  has  been  made.  The  original  speeific  aim  and  the  relevant  proposed 
work  for  eaeh  is  listed  below  and  addressed. 

Specifie  Aim  #1  stated:  “To  expand  and  refine  the  current  MIE  technique  to  enhance  its 
efficiency,  as  well  as  add  new  capabilities  such  as  handling  a  full  3D  or  combined  2D/3D 
elastodynamic  model  for  improved  accuracy.” 

An  improved  framework  is  in  progress  utilizing  parallel  processing  techniques.  In  order  to 
accommodate  the  methodology  of  MIE  in  creating  a  Jaeobian  matrix  fiilly  sensitive  to  the 
discretization  of  the  domain,  a  large  number  of  solutions  involving  the  finite  element  model  and 
the  subsequent  imge  deformation  are  required.  With  the  proposed  inerease  in  dimensionality,  the 
implementation  complexity  quiekly  increases  beyond  the  eapabilities  of  the  original 
MATLAB/FORTRAN/LAPACK  design.  Therefore,  the  Portable  Extensible  Toolkit  for 
Seientifie  Computation  (PETSc)  toolkit  [3,4]  was  seleeted  to  provide  the  neeessary  framework 
for  developing  sparse  matrix  system  solvers  and  split  the  Jacobian  formation  process.  A  separate 
C/C++  routine  has  also  been  written  to  perform  a  Gauss-Newton  optimization  and  interfaee  with 
PETSe  solver  structures.  A  current  typical  iteration  involves  the  solution  of  a  matrix  system  of 
approximately  le5  equations,  repeated  some  3000  times.  Utilizing  parallel  design  and  a  share  of 
100  CPUs  on  the  Vanderbilt  University  Advanced  Computing  Center  for  Researeh  and 
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Education  cluster,  this  has  been  tested  to  be  aehieved  on  the  order  of  30  minutes,  as  opposed  to 
original  estimates  on  available  sequential  processing  maehines  upwards  of  5000  minutes.  We 
note,  however,  that  to  effeetively  traverse  the  full  multi-dimensional  objeetive  function  space 
requires  several  (perhaps  tens  of)  iterations,  which  underscores  the  high  computational  demand 
of  the  method. 

We  have  preliminary  work  that  demonstrates  the  new  MIE  framework  in  action,  using  a  test  case 
of  simulated  deformation  based  on  elinieal  data.  Figure  1  shows  orthogonal  euts  in  the  three 
cardinal  anatomic  planes  for  an  image  volume  obtained  from  a  CT  scan  of  a  human  breast. 
Fibroglandular  tissue  can  be  visually  inspected  to  provide  eontrast  and  structure  from  adipose 
and  other  tissue  types.  The  test  case  involves  the  simulated  implantation  of  a  2-em  spherical 
tumor  at  the  eenter  of  the  breast  that  is  not  visible  within  the  intensity  field  of  the  image.  Guided 
by  a  finite  element  mesh  deformation  using  prescribed  boundary  displacements  (designed  to 
mimic  a  compression  source  as  deseribed  in  Speeifie  Aim  3),  a  target  image  volume  is  created. 
Diseretizations  of  the  model  and  image  domains  are  then  used  in  the  optimization  loop  to 
reeonstruct  the  inelusion.  Figure  2  shows  surfaee  renderings  of  the  image  volumes  and  the  results 
of  a  reeonstruetion  based  on  3200  spatially  distributed  elasticity  regions.  The  areas  delineated  to 
be  the  true  tumor  extent  contain  entites  that  the  algorithm  designated  as  having  stiffer  properties 
(~2x).  It  may  be  initially  seen  to  be  an  improper  chaeterization  given  that  the  faux  tumor  was 
actually  six  times  stiffer  than  the  surrouding  tissues.  However,  in  this  test  case,  the  inexact 
partitioning  of  the  mesh  elements  actually  caused  the  algorithm  to  search  for  a  tumor  about  3  cm 
in  size,  leading  to  a  compensatory  decrease  in  elastieity  contrast. 


Figure  1.  Orthogonal  cuts  of  a  CT  breast  scan  used  as  the  source  image  volume  for  simulation  studies.  From  left  to 
right:  axial,  coronal,  and  sagittal  views  as  designated  by  the  standard  anatomical  planes. 

In  order  to  further  explain  the  relationship  between  discretization  and  reeonstrueted  eontrast,  a 
partitioning  of  the  mesh  elements  utilizing  a  priori  knowledge  of  the  loeation  of  the  tumor  was 
used  to  elassify  the  tissue  types,  thereby  redueing  the  reeonstruetion  to  two  materials  instead  of 
the  3200.  This  reeonstruetion  very  favorably  followed  an  objective  funetion  minimization  to 
obtain  an  elastieity  contrast  of  the  inelusion  being  almost  exactly  six  times  stiffer  as  prescribed 
(see  Figure  3  below).  It  is  our  eurrent  assertion  that  shaping  the  objective  funetion  by 
dynamieally  rearranging  the  spatial  diseretization  of  the  model  during  the  optimization  ean  lead 
to  improved  elasticity  contrast  resolution,  and  studies  are  underway  to  address  this  issue. 
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Figure  2.  Top  row:  surface  renderings  in  perspective  of  the  breast  image  volume  in  both  the  native  undeformed 
(left)  and  deformed  (right)  states.  Bottom  row:  reconstruction  mappings  from  a  simulation  experiment.  Areas  of  high 
intensity  reflect  higher  (stiffer)  elasticity  values;  the  boundary  of  the  simulated  tumor  is  overlaid  by  the  circle. 


Proposed  work  involving  the  use  of  a  combined  2D-3D  model  is  under  investigation  but  not 
completed  at  this  time.  It  is  hoped  that  reducing  the  scope  of  the  problem  by  dimensionality  will 
facilitate  further  analysis  of  the  problem  by  alleviating  the  time  and  resource  restrictions  of  the 
full  implementation. 

Specific  Aim  #2  stated:  “To  perform  texture  analysis  on  input  images  in  order  to  quantify  a 
statistical  parameter  capable  of  estimating  the  success  of  elastographic  reconstruction.” 

Texture  analysis  and  noise  tolerance  testing  has  been  performed  with  statistical  quantification  of 
reconstruction  success.  Our  observations  during  the  ongoing  development  and  testing  of  the  MIE 
method  prompted  questions  concerning  the  quality  of  data  necessary  and  sufficient  to  achieve 
satisfactory  results  (i.e.,  the  fidelity  of  the  reconstructed  elasticity  image).  The  primary  inputs  to 
the  reconstruction  method  are  the  acquired  images  and  the  delineated  boundary  conditions  on  the 
domain.  While  it  is  clearly  preferable  to  have  idealized  data,  in  reality,  both  inputs  involve 
varying  levels  of  manual  interaction.  As  an  initial  study,  we  conducted  tests  on  the  effects  of 
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Figure  3.  Reconstruction  simulation  experiment  constrained  to  a  two-material  system  demonstrating  the  importance 
of  a  priori  information.  Knowing  the  location  of  the  inclusion  allowed  the  algorithm  to  quickly  search  the  objective 
function  space  and  arrive  (a)  back  at  the  original  [correct]  elasticity  contrast  (b).  Note  that  the  optimization  had 
actually  already  converged  at  iteration  7. 
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degradation  in  data  quality  on  the  end  reconstruction.  The  first  experiment  used  additive  image 
noise  to  obscure  the  underlying  texture  to  reflect  possible  scenarios  of  corruption  during 
acquisition.  Noise  fields  were  created  from  a  zero-mean  Gaussian  random  distribution  along  the 
variance  of  non-background  pixels  and  scaled  according  to  the  total  power  at  1,5,  10,  15,  20,  25, 
and  30%.  The  increasing  noise  had  a  confounding  impact  on  the  ability  of  the  similarity  metric  to 
make  a  proper  match.  It  was  found  that  the  reconstruction  was  tolerant  of  image  noise  up  to 
approximately  10%.  Figure  4  demonstrates  the  degradative  effects  of  image  noise. 


Figure  4.  (a)  Example  2D  reconstructions  resulting  from  the  distortion  of  the  target  image  using  additive  Gaussian 
random  noise  (from  top  left:  1,  5,  10,  20,  25,  30%).  The  true  elasticity  distribution  is  a  centrally  located  and  roughly 
circular  region,  and  the  noise  progressively  confounds  the  reconstruction,  (b)  The  decreasing  trend  of  reconstruction 
fidelity  as  determined  by  quantitative  evaluation  of  localizing  and  characterizing  the  detected  inclusion  from  a  given 
trial,  [see  Appendix  B] 

The  second  experiment  involved  boundary  condition  selection  error.  Currently,  point 
correspondence  at  the  outer  boundary  of  the  domain  is  determined  using  a  semi-automated  fit;  a 
polynomial  interpolation  is  used  to  make  the  initial  match,  and  then  corrected  manually  based  on 
salient  features.  By  perturbing  the  displacements  for  each  boundary  node  of  the  finite  element 
mesh,  typical  user  interaction  (e.g.  visual  cues  and  input  device  control)  can  be  simulated.  A  gold 
standard  set  of  boundary  conditions  known  to  produce  an  accurate  reconstruction  was  modified 
using  randomized  vectors  of  equal  magnitude  (0.1,  0.2,  0.3,  0.5,  0.75,  1.0,  1.5,  and  2.0  pixel 
units)  reflecting  a  range  of  typical  localization  skill  for  users  from  poor  to  expert.  It  was  noted 
that  randomizing  all  vectors  can  actually  result  in  twisting  of  elements  that  resulted  in  significant 
alterations  of  displacements  in  the  interior  of  the  mesh,  leading  to  grossly  inaccurate  model 
deformations.  However,  boundary  condition  mismapping  of  less  than  0.5  pixel  units  was 
generally  tolerated  by  the  reconstruction  algorithm,  a  range  that  is  relatively  easily  achieved  for 
most  users.  For  further  detail,  please  see  the  full  text  of  this  work  as  listed  in  Appendix  B. 

It  should  be  noted  that  the  ordinate  axis  in  Figure  4(b)  is  a  ‘quality  of  reconstruction  score’ 

(QRS)  that  has  been  developed  within  the  context  of  this  project  in  order  to  quantify  the 
localization  accuracy  of  the  method.  In  comparison  and  conjuction  with  more  standard  measure 
in  the  elastography  field  of  contrast-to-noise  ratio  [5],  QRS  has  been  used  to  determine  relevant 
positional  and  material  characterization  in  both  simulation  and  data  studes.  The  metric  is 
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determined  by  a  classification  of  the  reconstruction  [6]  that  is  then  compared  to  the  (known) 
segmentation  of  the  actual  elasticity  distribution.  By  examining  the  rate  of  accurately  selecting 
an  element  of  the  domain  to  be  of  the  correct  material  class,  a  conditional  probability  closely 
related  to  the  positive  predictive  value  of  the  test  is  obtained;  we  have  determined  a  posteriori 
that  a  QRS>80%  is  typically  indicative  of  a  successful  reconstruction.  The  use  of  QRS  can  (and 
will  be)  similarly  applied  to  the  analysis  of  forthcoming  3D  reconstuction  experiments.  For 
further  discussion,  please  see  material  in  Appendix  A  and  B. 

Proposed  work  involving  the  use  of  a  feature  tracking  and  frequency  domain  analysis  is  under 
investigation  but  not  completed  at  this  time.  As  more  data  sets  are  collected,  it  is  hoped  that 
establishing  a  pattern  for  understanding  the  reconstruction  algorithm  behavior  will  become 
statistically  relevant. 

Specific  Aim  #3  stated:  “To  engineer  a  device  that  can  accurately  produce  compressive 
forces  necessary  for  phantom  setups  within  current  clinical  imaging  systems,  providing  the 
basis  for  a  future  device  that  can  be  used  in  a  clinical  setting.” 

A  compression  device  has  been  constructed  and  tested  in  magnetic  resonance  (MR)  and  X-ray 
computed  tomography  (CT)  imaging  systems  using  a  polyvinyl  alcohol  phantom  and  contrast 
agents.  The  compression  device  is  composed  of  a  rectangular  Plexiglas  frame  that  traps  the 
phantom  in  at  least  two  directions  with  a  sliding  wall  and  the  compression  plate,  which  houses  an 
air  bladder  in  a  polycarbonate  frame.  When  inflated,  the  air  bladder  provides  a  deformation  of  up 
to  5  cm.  The  prototypical  phantom  used  has  been  fabricated  as  a  polyvinvyl  alcohol  cryogel 
(~650  cc,  6-8%  wt/vol,  1  or  2  free-thaw-cycles)  in  a  manner  consistent  with  the  methods 
presented  in  [7,  8].  The  result  is  a  dome  shape  approximately  10-11  cm  at  the  base  and  tapering 
over  a  depth  of  about  5-6  cm.  The  system  has  been  imaged  by  both  MR  and  CT  scans,  with 
contrast  agents  of  copper  sufate  and  iodine  solution,  respectively,  having  been  used  to  enhance 
the  signals.  Volumetric  renderings  of  example  scans  using  minimal  post-processing  are  shown 
below  in  Figure  5;  note  that  the  device  frame  is  clearly  present  in  the  CT  image  because  of  its 
density  but  is  invisible  on  MR. 

A  prototype  compression  chamber  that  is  more  clinically  oriented  has  been  designed  to  fit  into 
the  chassis  of  a  Philips  Intera  breast  coil  unit.  It  has  just  recently  been  fabricated  from  clear 
acrylic  tube  segments  in  which  the  air  bladders  are  attached  using  polycarbonate  pins  and  then 
covered  with  an  expandable  nylon  sheet  (see  Figure  6).  Image  acquisition  studies  are  currently 
being  designed  to  utilize  this  system. 
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Figure  5.  Top  row:  Photographs  of  the  polyvinyl  alcohol  phantom  inside  the  compression  device  without  (left)  and 
with  compression  (right).  Bottom  row:  Surface  renderings  of  from  image  volumes  obtained  on  a  breast  phantom  in 
MR  (left)  and  CT  (right)  scanners  while  enclosed  in  the  compression  device  and  with  the  air  bladder  engaged  to 
provide  a  deformation. 


Figure  6.  Left:  Schematic  of  compression  device  designed  for  clinical  use  breast  coil.  Right:  Photograph  of 
assembly  looking  down  into  the  Philips  Intera  chassis. 
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Key  Research  Accomplishments 


It  has  been  observed  both  prior  to  and  during  the  attempt  to  extend  the  MIE  technique  to  cover 
fully  three-dimensional  data  that  the  problem  is  potentially  intractable  due  to  the  large 
computational  demand  and  inherently  ill-constrained  system.  The  demonstration  of  a  parallelized 
code  base  is  for  this  project  a  significant  finding  in  that  it  confirms  that  operating  on  volumetric 
images  and  models  is  a  reasonable  working  hypothesis,  albeit  still  challenging. 


Reportable  Outcomes 

Work  on  the  MIE  method  has  so  far  resulted  in  two  conference  papers  and  an  additional  poster 
presentation.  Prior  work  that  was  completed  in  the  reporting  timeline  resulted  in  a  peer-reviewed 
journal  publication.  These  items  are  in  part  providing  the  foundation  for  a  thesis  proposal  to  be 
submitted  in  the  near  future.  Didactic  coursework  requirements  for  the  PhD  degree  have  also 
been  completed  at  this  time. 

Poster  presentations 

Vanderbilt  University  Institute  of  Imaging  Science  retreat  (June  2005) 

Vanderbilt  University  Medical  Scientist  Training  Program  retreat  (July  2005) 

Conference  papers 

Ou  JJ,  Barnes  SL,  Miga  MI,  “Application  of  multi-resolution  modality  independent 
elastography  for  detection  of  multiple  anomalous  objects”,  SPIE  Medical  Imaging  2006. 

Ou  JJ,  Barnes  SL,  Miga  MI,  “Preliminary  testing  of  sensitivity  to  input  data  quality  in  an 
elastographic  reconstruction  method”,  IEEE  International  Symposium  on  Biomedical  Imaging 
2006. 

Schuler  DR,  Ou  JJ,  Barnes  SL,  Miga  MI,  “Automatic  surface  correspondence  methods  for  a 
deformed  breast”,  SPIE  Medical  Imaging  2006. 

Journal  publications 
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Conclusions 

The  current  results  and  progress  denoted  in  this  report  are  within  the  proposed  statement  of  work 
and  are  encouraging  towards  completion  of  the  overall  objectives  with  further  effort.  No 
significant  deviations  are  reported  at  this  time. 
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ABSTRACT 

This  work  extends  a  recently  realized  inverse  problem  technique  of  extracting  soft  tissue  elasticity  information  via  non- 
rigid  model-based  image  registration.  The  algorithm  uses  the  elastic  properties  of  the  tissue  in  a  biomechanical  model  to 
achieve  maximal  similarity  between  image  data  acquired  under  different  states  of  loading.  A  new  multi-resolution,  non¬ 
linear  optimization  framework  has  been  employed  which  allows  for  improved  performance  and  object  detection.  Prior 
studies  have  demonstrated  successful  reconstructions  from  images  of  a  tissue-like  thin  membrane  phantom  with  a  single 
embedded  inclusion  that  was  significantly  stiffer  than  its  surroundings.  For  this  investigation,  a  similar  phantom  was 
fabricated  with  two  stiff  inclusions  to  test  the  effectiveness  of  this  method  in  discriminating  multiple  smaller  objects. 
Elasticity  values  generated  from  both  simulation  and  real  data  testing  scenarios  provided  sufficient  contrast  for  detection 
and  good  quantitative  localization  of  the  inclusion  areas. 

Keywords:  Elastography,  elasticity  imaging,  multi-resolution  methods,  image  similarity,  finite  elements 

1.  INTRODUCTION 

The  practice  of  palpating  soft  tissue  structures  in  the  course  of  the  clinical  physical  exam  has  had  a  long-standing 
history  of  providing  correlation  of  improper  stiffness  with  pathology.  The  ability  to  characterize  the  mechanical 
properties  of  tissue  is  a  potential  source  of  additional  information  relevant  for  detection  and  diagnosis  of  a  disease 
process,  and  has  implications  for  the  assessment  of  treatment.  One  way  in  which  this  could  be  achieved  in  a  minimally 
invasive  manner  is  by  analyzing  tissue  deformation  through  imaging  and/or  image  processing  techniques,  which  is  a 
central  goal  of  the  field  of  elastography  [1].  Application  of  such  methods  to  the  interrogation  of  the  breast  [2,3],  skin 
[4-6],  prostate  [7],  and  other  accessible  organ  systems  is  an  emerging  area  of  research. 

Many  of  the  current  elastography  methods  are  founded  in  ultrasound  (US)  and  magnetic  resonance  imaging 
(MR)  and  involve  the  estimation  of  induced  displacements  within  the  tissue  of  interest  to  infer  the  elasticity  distribution. 
We  have  pursued  the  development  of  a  reconstruction  method  utilizing  quasi-static  deformation  and  image  similarity 
metrics  that  has  been  termed  'modality-independent  elastography’  (MIE)  [8-10]  because  of  its  potential  to  handle  native 
anatomical  image  data  from  different  modalities  with  simple  modification  to  the  acquisition  procedure.  Common 
problems  facing  all  of  these  methods  involve  limitations  with  the  accurate  recovery  of  elastic  property  values,  detection 
of  small  lesions  in  tissue,  and  the  resolution  of  multiple  discrete  lesions  [11,12].  Building  upon  recent  study  involving  a 
single  focal  lesion  [6],  the  objectives  of  this  work  were  to  challenge  the  ability  of  the  MIE  method  to  reconstruct  a 
scenario  of  two  small  inclusions  embedded  in  a  homogeneous  domain  and  to  further  explore  the  feasibility  of  the 
method  in  handling  image  data  from  different  imaging  modalities.  This  was  accomplished  by  performing  simulated 
reconstructions  using  images  obtained  from  X-ray  computed  tomography  (CT),  MR,  and  digital  photography  and  then  a 
reconstruction  from  a  real-world  experiment  using  a  thin  phantom  membrane. 

2.  METHODS 


2.1  Elastographic  reconstruction  framework 

The  conceptual  framework  for  our  elastographic  reconstruction  has  been  previously  described  in  [6,8-10].  In  brief,  an 
image  of  a  tissue  of  interest  {source)  is  deformed  by  a  biomechanical  computer  model  and  compared  against  an 
acquired  image  of  the  same  tissue  in  a  mechanically  loaded  state  {target).  The  deformation  and  comparison  is  repeated 
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using  systemic  updates  of  elasticity  parameters  until  a  suitable  match  in  intramodal  image  similarity  is  achieved  in  a 
least  squares  manner  to  satisfy  a  multi-resolution,  non-linear  optimization  scheme.  This  process  can  be  classified  as  an 
inverse  problem,  with  model-based  deformation  of  the  source  image  representing  the  forward  problem.  Each  of  the 
three  major  components  (model,  image  comparison,  and  optimization)  is  described  in  more  detail  in  the  following 
sections,  and  a  flow  chart  representation  of  the  overall  process  is  included  in  Figure  1. 

2.1.1  Biomechanical  model 

A  central  component  to  the  model-based  inverse  problem  is  the  manner  in  which  the  continuum  is  represented.  While 
the  constitutive  model  that  best  describes  tissue  deformation  mechanics  is  more  complex,  for  this  work,  linear  isotropic 
elasticity  has  been  employed.  The  partial  differential  equation  that  expresses  a  state  of  mechanical  equilibrium  can  be 
written  as  [13]: 

V^a  =  0  (1) 


where  cris  the  Cartesian  stress  tensor. 

For  the  purposes  of  the  following  experimentation,  we  also  apply  either  the  plane  stress  or  plane  strain 
approximations  to  the  thin  membrane  and  breast  cross-section  trials,  respectively.  The  direct  consequence  of  this  is  a 
reduction  of  the  36  stiffness  constraints  in  the  general  3D  formulation  of  Cauchy’s  Faw  to  the  two  parameters  of 
Young’s  modulus  (E)  and  Poisson’s  ratio  (v)  in  2D.  These  simplifications,  while  significant,  are  appropriate 
descriptions  of  sufficiently  thin  and  thick  systems  under  planar  loading.  In  plane  stress. 
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describes  the  constitutive  relationship  between  the  Cartesian  stress  tensor  [a^,  Qy,  Xxy]  and  strain  tensor  [Sx,  Sy,  yxy]. 
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A  finite  element  (FE)  model  using  triangular  elements  is  constructed  from  the  source  image  and  assigned  appropriate 
boundary  conditions  based  on  estimated  displacement  or  stress  (i.e.  Dirichlet  and  Neumann  conditions,  respectively). 
The  standard  Galerkin  method  of  weighted  residuals  [14]  is  used  to  construct  and  solve  the  system. 

2.1.2  Image  deformation  and  comparison 

To  further  describe  the  reconstruction  process,  we  introduce  some  additional  terminology  at  this  point.  The  model 
domain  is  equivalent  to  the  total  area  of  the  FE  mesh  constructed  using  the  source  image  as  stated  above  and  contains 
the  relevant  elasticity  information.  The  model  domain  is  partitioned  by  a  K-means  clustering  of  the  element  centroids 
(MATFAB  R14,  Mathworks,  Natick,  MA)  into  N  number  of  regions,  each  of  which  has  a  distinct  set  of  spatially 
homogeneous  elastic  properties.  Subdividing  in  this  manner  allows  for  the  implementation  of  the  multi-resolution 
reconstruction  whereby  progressively  finer  spatial  distributions  of  elasticity  parameters  are  utilized  in  the  process,  a 
method  that  improves  upon  previous  versions  using  only  a  single  resolution  [8-10].  Analogously,  the  comparison 
domain  is  an  area  specified  by  semi-automated  segmentation  on  the  target  image  and  contains  information  pertaining  to 
image  similarity.  The  comparison  domain  is  separated  into  M  number  of  rectangular  zones  containing  approximately 
equal  numbers  of  pixels. 
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The  reconstruction  algorithm  begins  by  assigning  an  initial  Young’s  modulus  value  to  each  of  the  regions  at 
the  coarsest  resolution.  Poisson’s  ratio  is  held  constant  at  v  =  0.485  to  represent  a  nearly  incompressible  material.  The 
FE  model  is  solved  to  determine  the  nodal  mesh  displacements,  which  are  in  turn  used  to  deform  the  source  image.  This 
model-deformed  image  is  then  compared  to  the  target  image  for  every  zone  using  an  intensity-based  image  similarity 
metric.  While  a  number  of  methods  are  available  for  such  a  task,  here,  we  utilize  the  correlation  coefficient  (CC)  [15] 
throughout,  as  it  has  empirically  demonstrated  superior  performance  over  other  metrics  such  as  the  sum  of  squared 
differences  and  normalized  mutual  information. 


2.1.3  Optimization  scheme 


Let  r  be  a  function  that  represents  the  model-based  image  deformation  and  takes  as  its  input  a  vector  of  elastic  modulus 
values  E  of  length  N  that  corresponds  to  the  current  distribution  of  regions  in  the  model  domain.  Then  for  two 
distributions  of  modulus  values  Ei  and  E2,  the  similarity  between  the  images  produced  by  r(Ei)  and  r(E2)  is  the  vector 
S  of  length  M  containing  evaluations  of  the  correlation  coefficient  corresponding  to  the  distribution  of  zones  in  the 
comparison  domain.  The  elasticity  parameter  optimization  can  be  written  as  the  minimization  of  the  least  squares  error 
objective  function 


W  —  K  ^ 

^  \^TRUE  ^EST\ 
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where  Strue  is  the  set  of  similarity  values  achieved  when  comparing  the  target  image  to  itself,  Srst  is  the  similarity 
between  the  model-deformed  source  and  the  target  images  using  current  estimates  of  the  elastic  modulus  distribution, 
and  !•!  denotes  the  vector  L2  norm.  By  definition,  Strue  is  the  maximum  value  for  the  similarity  metric  (max  CC  =1). 
Using  a  Levenberg-Marquardt  approach,  the  residual  form  of  equation  (4)  becomes 

[jV  +  a/]{A£}  =  [/T'S™j 

where  J  =  ^Sest/^E  is  the  Jacobian  matrix  of  size  MxN  and  I  is  the  TV x  identity  matrix.  Because  J^J  is  typically  an 
ill-conditioned  term,  the  regularization  parameter  a  is  determined  using  the  methods  described  in  [16].  Modulus  values 
of  the  regions  at  a  given  resolution  are  updated  by  AE  until  an  error  tolerance  is  reached  or  a  maximum  number  of 
iterations  have  been  completed.  Upon  reaching  a  stopping  criterion,  the  material  property  description  is  interpolated 
onto  the  next  (i.e.  finer)  resolution  and  the  above  steps  are  repeated.  Spatial  averaging  of  modulus  values  within  the 
model  domain  and  solution  relaxation  between  successive  resolution  levels  are  also  utilized  to  improve  the  stability  of 
the  optimization. 
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Figure  1.  Flow  chart  of  elastographic  reconstruction  framework. 

2.2  Reconstruction  experiments 

A  two-material  phantom  membrane  of  simulated  skin  had  been  previously  constructed  [6]  using  Smooth-On™ 
polyurethanes  (Smooth-On,  Easton,  PA)  designated  by  the  manufacturer  as  Evergreen  10  and  Evergreen  50.  These 
materials  have  essentially  indistinguishable  colors  but  vary  significantly  in  their  elastic  modulus  values,  so  the  former 
was  used  as  the  bulk  material  and  the  latter  for  stiff  objects.  From  material  testing,  the  elastic  modulus  contrast  was 
expected  to  be  approximately  5.7:1.  The  phantom  was  made  to  contain  two  circular  stiff  inclusions  1.5  cm  in  diameter 
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embedded  near  opposing  comers  of  a  rectangular  field  of  bulk  material  measuring  15  cm  x  14  cm.  A  black  permanent 
marker  was  used  to  place  a  pattern  of  regularly  spaced  (~1  cm)  grid  lines  across  the  membrane.  The  thin  membrane  was 
securely  clamped  along  two  opposite  edges  and  then  subjected  to  a  uniaxial  tensile  displacement  (-8%  strain)  by  means 
of  a  milling  vise.  A  commercial  webcam  (Logitech  QuickCam  Pro  4000,  960  x  1280  pixel  resolution)  was  rigidly 
mounted  above  the  membrane  to  acquire  image  pairs  of  the  pre-  and  post-stretched  states. 

To  initially  test  the  method  regarding  the  two-inclusion  scenario,  a  simulation  using  the  source  image  of  the 
membrane  was  performed  by  deforming  it  with  a  prescribed  model  (plane  stress)  of  known  boundary  displacements  and 
elasticity  parameters  to  generate  a  target  image;  high  modulus  values  were  assigned  to  elements  bounded  by  a 
segmentation  of  the  inclusion  locations.  A  reconstruction  was  then  performed  using  the  actual  image  data  acquired  as 
described  above.  In  both  cases,  resolutions  of  A=  16,  64,  256,  512,  and  800  regions  and  M=  400  zones  were  used.  The 
results  of  the  idealized  and  real  data  reconstructions  are  shown  in  Figures  4  and  5,  with  further  quantitative  evaluation  in 
Table  1. 

In  order  to  examine  the  robustness  of  the  method  regarding  its  use  of  data  from  differing  sources,  simulation 
reconstructions  were  performed  using  image  slices  extracted  from  breast  image  volumes  obtained  from  CT  and  MR 
scans  (see  Figure  3).  Although  these  were  taken  from  two  different  patients,  the  images  were  selected  to  be 
approximately  corresponding  slices  ~2  cm  away  from  the  chest  wall  in  the  coronal  orientation  of  the  standard 
anatomical  position.  The  simulations  were  set  up  in  the  same  manner  as  for  the  digital  photographs,  using  either  one  or 
two  inclusions  of  about  1  cm  in  diameter  embedded  within  the  true  elasticity  distribution  and  a  small  compression  (~8% 
strain)  in  the  cranial-caudal  direction.  The  relative  stiffness  of  the  inclusions  was  designated  to  be  5.7:1  for  consistency 
with  the  material  testing  data  and  also  because  the  value  is  fairly  representative  of  breast  tumor  properties  [17].  The 
plane  strain  model  approximation  was  used  in  the  breast  simulation  trials,  progressing  through  resolutions  of  A  =  24,  64, 
256,  and  576  regions  using  M  =  200  zones.  The  reconstruction  method  was  then  run  for  all  four  test  cases,  and  the 
results  are  presented  in  Figures  6  and  7  and  Table  2. 


Figure  2.  (Left  to  right):  Phantom  membrane  in  undeformed  state  {source  image),  under  deformation  {target  image),  and  differenee 
image.  Arrows  in  the  left  panel  indieate  the  positions  of  the  two  stiff  inelusions. 


2.3  Reconstruction  evaluation 

The  fidelity  of  the  elasticity  reconstruction  was  evaluated  on  its  ability  to  detect  the  presence  of  an  inclusion  based  on 
classification  of  the  material  property  distribution,  and  the  retrospective  accuracy  of  localizing  the  lesions.  The  elastic 
properties  as  a  whole  were  treated  as  a  Gaussian  mixture  of  two  classes  and  separated  by  a  threshold  established  via  the 
method  described  in  [18].  The  likelihood  of  detecting  a  lesion  in  the  elasticity  image  was  found  using  the  contrast-to- 
noise  ratio  as  defined  by  [12,19]: 
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Figure  3.  Images  sliees  of  breast  tissue  extraeted  from  a  CT  volume  (left)  and  MR  volume  (right)  used  in  simulation  study  of  the 
ability  of  the  reeonstruetion  method  to  utilize  disparate  image  data  types. 


where  //  and  <j^  are  the  sample  mean  and  variance  of  a  material  property  distribution  and  the  subscripts  L  and  B  denote 
the  lesion  and  bulk  material  classes,  respectively.  As  a  quantitative  assessment  of  the  localization  of  the  lesion(s),  the 
positive  predictive  value  of  correctly  identifying  a  lesion  material  within  the  known  segmented  region  of  the  inclusions 
was  used  as  a  'quality  of  reconstruction  score'  (QRS).  This  value  is  significant  because  identification  of  the  lesion 
border  and  material  classification  are  done  independently,  so  any  user  knowledge  of  the  test  scenario  does  not  influence 
the  performance  of  the  measure.  Cutoffs  for  successful  detection  and  localization  were  set  at  CNR>2.2  as  noted  by  [12] 
and  QRS>80%  as  determined  by  prior  study  in  our  laboratory.  The  average  modulus  contrast  is  found  from  the  ratio  of 
the  means  of  the  two  material  classes,  and  a  peak  modulus  contrast  value  is  also  reported  by  taking  the  ratio  of  two 
manually  selected  homogeneous  regions  of  approximately  equal  area  known  to  be  representative  of  the  two  materials. 
It  should  be  noted  that  in  other  work  not  presented  here,  the  definition  of  QRS  included  a  weighting  factor  provided  by 
the  estimated  reconstruction  modulus  contrast,  but  for  the  current  purposes,  only  localization  accuracy  was  considered 
to  maintain  an  objective  evaluation  of  inclusion  detection. 


3.  RESULTS 

Figure  4  demonstrates  the  ability  of  the  reconstruction  method  to  produce  an  elasticity  map  from  the  simulation  data 
with  good  localization  of  the  inclusions  that  are  easily  visually  distinguishable  from  the  surrounding  bulk  material.  The 
progression  through  resolutions  of  N  =  64,  256,  512,  and  800  regions  shows  improving  delineation  of  the  inclusions  and 
elastic  contrast.  Figure  5  demonstrates  a  similar  behavior  for  the  reconstruction  of  the  acquired  phantom  membrane 
data,  with  both  spatial  definition  and  modulus  contrast  increasing  with  the  finer  discretization.  Table  1  summarizes  the 
quantitative  evaluation  of  the  reconstructions  in  both  simulation  and  phantom  trials,  including  CNR,  contrast  ratio,  and 
QRS  values.  The  CNR  values  are  sufficient  to  allow  for  discrimination  of  the  two  materials  and  the  identification  of  the 
inclusions  was  determined  to  be  accurate  in  both  cases.  The  reconstruction  of  the  phantom  membrane  does  show  some 
misclassification  along  the  border  where  the  deformation  was  applied  as  well  as  in  the  comer  adjacent  to  one  of  the 
inclusions  (see  Figure  5d). 

Figures  6  and  7  show  the  final  reconstmction  results  for  the  CT  and  MR  breast  slice  simulations  using  either 
one  or  two  inclusions.  In  both  test  scenarios,  the  resolvability  of  the  stiffer  material  was  found  to  be  adequate  according 
to  the  CNR  threshold,  but  definitely  higher  in  the  MR-derived  elasticity  images.  Localization  of  the  inclusions  yielded 
excellent  QRS  values  in  reconstmctions  using  either  modality,  again  higher  (though  slightly)  for  the  MR  images. 
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Figure  4.  Reconstruction  of  the  simulated  membrane  deformation  using  idealized  model  parameters,  progressing  through  finer 
resolution  distributions  (a)-(d)  of  64,  256,  512,  and  800  regions. 


Figure  5.  Reconstruction  of  the  actual  membrane  data.  A  faint  contour  in  (d)  is  present  to  demarcate  the  actual  position  of  the  stiff 
inclusions.  Again,  panels  (a)-(d)  demonstrate  the  effect  of  the  multi-resolution  method  in  utilizing  64,  256,  512,  and  800  regions  to 
better  capture  the  shape  and  location  of  the  inclusions. 


Table  1.  Quantitative  reconstruction  evaluations. 


Avg  CR 

Max  CR 

CNR 

QRS  (%) 

Simulation 

2.7 

4.0 

4.4 

97.7 

Phantom 

2.6 

4.1 

2.8 

88.5 
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Figure  6.  Reconstructions  of  simulation  trials  for  the  CT  breast  slice  using  a  single  inclusion  (left)  and  two  inclusions  (right).  The 
true  inclusion  boundaries  are  overlaid  in  each  elasticity  image. 


Figure  7.  Reconstructions  (bottom  row)  of  simulation  trials  for  the  MR  breast  slice  using  a  single  inclusion  (left)  and  two  inclusions 
(right).  The  true  elasticity  distributions  are  also  shown  (top  row)  for  comparison. 


Table  2.  Quantitative  reconstruction  evaluations. 


Avg  CR 

MaxCR 

CNR 

QRS  (%) 

CT  (1  inclusion) 

2.1 

3.1 

3.0 

97.6 

CT  (2  inclusions) 

2.0 

2.6 

3.5 

96.9 

MR  (1  inclusion) 

2.8 

3.7 

20.0 

100 

MR  (2  inclusions) 

2.7 

3.7 

5.7 

99.8 
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4.  DISCUSSION 

The  results  of  the  phantom  membrane  experiment  are  encouraging  because  of  their  similarity  to  the  idealized 
simulation.  Despite  nonlinear  model-data  mismatch,  out-of-plane  distortions  during  stretching,  and  possible  boundary 
condition  inaccuracies,  the  elasticity  reconstruction  demonstrated  good  localization  of  the  two  small  inclusions.  The 
majority  of  the  problems  in  reconstruction  are  mostly  likely  due  to  noise  incurred  in  the  mapping  of  the  boundary 
displacements.  It  should  be  noted  that  the  phantom  reconstruction  was  achieved  with  a  non-pigmented  lesion  (see 
Figure  2,  arrows),  indicating  that  deflections  of  the  image  structure  are  capable  of  driving  the  image  similarity  metric  of 
the  reconstruction  process.  This  does  intuitively  suggest  that  some  metric  for  rating  the  complexity  and  density  of 
image  pattern  in  relation  to  algorithm  success  may  be  important  and  is  currently  under  investigation.  Preliminary  data 
not  presented  in  this  work  indicates  that  such  a  threshold  does  exist  for  image  data  that  can  be  properly  analyzed  by  the 
current  framework.  The  modality  independence  of  the  method  is  also  supported  by  the  results  here;  clearly,  the 
Hounsfield  units  of  CT,  floating  point  values  from  an  MR  volume,  and  the  luminance  captured  by  the  CCD  sensor  of  a 
digital  camera  are  quite  different  types  of  data  to  handle  because  they  are  based  on  different  physical  principles.  The 
simulation  reconstructions  demonstrate  that  the  method  is  indifferent  to  these  differences  by  treating  the  data  as  an 
arbitrary  range  of  intensities  and  will  converge  towards  the  true  elasticity  distribution  based  on  the  image  pattern 
available.  This  is  a  possible  explanation  for  the  qualitatively  more  satisfactory  results  from  the  MR  simulations 
compared  to  the  CT  trials  because  the  distribution  of  intensities  from  the  former  modality  yielded  a  more  diversified 
histogram,  an  attribute  that  should  naturally  aid  an  intensity-based  metric. 

While  an  ideal  reconstruction  would  also  be  accurate  in  characterizing  a  lesion  by  its  modulus  contrast,  our 
focus  in  the  study  was  to  test  the  ability  of  the  method  to  detect  and  localize  the  inclusions.  In  previous  experimentation 
with  reconstructions  of  single  focal  lesions,  we  have  been  generally  successful  in  achieving  a  contrast  ratio  within  25% 
of  the  true/expected  value.  It  is  somewhat  troubling  that  the  contrast  ratios  calculated  here  did  not  meet  that  criterion, 
although  the  experiments  with  the  phantom  membrane  came  fairly  close  (28%).  However,  these  results  underscore  the 
difficulty  of  the  scenarios  in  not  only  having  to  deal  with  multiple  inclusions  but  quite  small  ones  in  both  the  true 
physical  sense  and  also  the  scale  of  the  domain.  Any  of  the  given  inclusions  tested  in  simulation  and  with  the  real  data 
were  detected  within  a  homogeneous  domain  approximately  an  order  of  magnitude  larger  (e.g.,  1.5-cm  lesions  in  a  15 
cm  X  14  cm  domain  for  the  phantom).  The  expectation  of  being  able  to  identify  with  any  confidence  the  presence  of  the 
inclusion  is  comparable  to  the  observations  made  in  [12]  where  the  test  of  finding  a  single  5 -mm  lesion  within  a  4  cm  x 
5  cm  domain  proved  to  be  the  most  problematic.  Therefore,  the  localization  of  the  lesions  as  determined  by  the  CNR 
and  QRS  metrics  is  deemed  to  be  a  success,  and  further  investigation  into  the  nature  of  the  method  with  respect  to  the 
scale  of  the  lesion  and  domain  is  warranted. 


5.  CONCLUSIONS 

In  this  work,  we  have  presented  further  testing  of  a  method  for  recovering  elasticity  parameters  by  maximizing  the 
similarity  between  images  of  a  tissue  of  interest  acquired  under  two  different  states  of  quasi-static  loading  within  the 
context  of  an  inverse  problem.  The  specific  experiments  presented  here  examined  the  effectiveness  of  the  technique  for 
the  detection  of  multiple  small  discrete  focal  lesions  embedded  in  an  otherwise  homogeneous  medium,  as  well  as 
further  proof-of-concept  work  in  its  applicability  to  utilize  image  data  from  various  modalities.  In  both  cases,  the 
method  provided  accurate  localization  of  the  lesions  based  on  the  reconstruction  of  relevant  elasticity  contrast.  Because 
the  biomechanical  model,  multi-resolution  optimization,  and  image  acquisition  are  each  modular  components  of  the 
framework,  this  elastographic  reconstruction  technique  is  readily  extensible  for  added  sophistication,  and  there  is 
ongoing  work  to  enhance  the  methodology  with  more  complex  models  and  advances  in  imaging  technology. 
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Modality  Independent  Elastography 
(MIE)  Concepts 


•  (Solid)  tumors  are  usually  stiffer  than  surrounding 
tissue 


•  Soft  tissue  interrogation  of  various  organ  systems 
(e.g.  skin,  liver,  prostate,  breast)  for  tumor 
detection 


•  Elastography  gives  representation  of  a  structure 
according  to  its  mechanical  properties 

•  Deformation  processes  indicative  of  material 
inhomogeneity  can  be  captured  by  imaging  and 
approximated  with  modeling 

•  Associate  form  and  function  through  image  analysis 
separate  from  modality  acquisition 
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MIE  Components 

•  (1)  Biomechanical  FE  model  of 

soft-tissue  deformation 

•  Conservation  of  stresses  (continuum) 

•  Constitutive  stress-strain  relation  (Hooke's  Law) 


<j  =  Es 


MIE  Components  (cont.) 

•  (2)  Similarity  measure  for 

comparing  images 

•  Acquired  "pre-"  (source)  &"post-"  (target) 
quasi-static  deformation 

•  Intensity-based  registration  metrics 
•  MI,  NMI,  SSD,  CC,  GC 
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MIE  Components  (cont.) 


•  (3a)  Optimization  routine  to  update 
material  properties  in  the  model 


•  Objective  function  based  on  similarity 


MIE  Components  (cont.) 

•  (3b)  Discretization  of  eiasticity 
distribution  and  image  data 

•  Multi-resolution  K-means  clustering  of  elements 
("regions") 

•  Sampling  of  image  comparison  area  ("zones") 
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MIE  Framework 


Medical  Physics,  vol.  32,  no.  5,  pp.  1308-1320,  2005 
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Study  Objectives: 

Further  Testing  of  MIE 

•  Modality  independence 

Digital  photography 

X-ray  computed  tomography  (CT) 

Magnetic  resonance  (MR) 

•  Two  (small)  inclusions 

•  Simulation  and  phantom  membrane 
study 
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Evaluating  MIE 

•  Classify  reconstruction 

Two-class  Gaussian  mixture  model 

•  Detectability  via  elasticity  image  contrast 


•  Localization  accuracy 

Positive  predictive  value  of  identifying  lesion 
material  in  correct  location 


CT  breast  slice  -  simulation 

Source  Targetl  Target2 
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MR  breast  slice  -  simulation 

Source  Targetl  Target2 
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Summary 

•  Modality  independence  via  simulation  for 
handling  various  data  types 

•  Multi-resolution  approach  potentially 
improves  optimization  convergence 

•  Two  small  stiff  inclusions  reconstructed  in 
phantom  membrane  experiment 

•  Detectability  accomplished  via  CNR 

•  Localization  successful  as  evaluated  by  QRS 
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<finis> 


Future  Directions  /  Research  Questions 

•  Biological  tissues  are  not  typically  linearly 
elastic 

•  Need  for  accurate  boundary  conditions 
creates  dependence  on  segmentation 
methodology 

•  Not  all  data  sets  necessarily  contain 
sufficient  information  for  elastographic 
reconstruction 
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ABSTRACT 

An  elastographic  reconstruction  method  has  been  developed 
to  recover  the  material  properties  of  soft  tissue  by  model- 
based  analysis  of  image  data  acquired  at  different  states  of 
mechanical  loading.  The  algorithm  utilizes  image  similarity 
as  part  of  the  cost  function  for  a  multi-resolution,  non-linear 
optimization.  Previous  work  with  a  phantom  membrane 
used  for  simulated  dermoscopic  application  has  prompted 
this  preliminary  investigation  of  the  relative  effects  of 
additive  image  noise  and  boundary  condition  determination 
errors  on  the  performance  of  the  method.  The  results  as 
quantified  by  elasticity  contrast  and  localization  accuracy 
indicate  that  the  reconstruction  process  is  robust  in  the 
presence  of  realistic  levels  of  image  corruption  and  tolerates 
the  majority  of  boundary  condition  mapping  errors. 


1.  INTRODUCTION 

The  practice  of  palpating  soft  tissue  structures  in  the  course 
of  the  physical  exam  for  assessing  tissue  health  has  had  a 
long-standing  clinical  history  of  providing  correlation 
between  improper  stiffiiess  and  pathology.  The  ability  to 
characterize  the  mechanical  properties  of  tissue  is  therefore 
a  potential  source  of  information  relevant  for  both  diagnosis 
and  prognosis.  One  way  in  which  this  could  be  achieved  in  a 
non-invasive  manner  is  through  analysis  of  tissue 
deformation  with  imaging  and  image  processing  techniques, 
which  is  a  central  goal  of  the  field  of  elastography  [1]. 

The  conceptual  framework  for  our  elastographic 
reconstruction  has  been  previously  described  in  [2-4].  In 
brief,  images  of  a  tissue  of  interest  are  acquired  in  an  initial 
{source)  and  then  mechanically  loaded  state  {target).  The 
source  image  is  deformed  by  a  prescribed  computational 
model  and  compared  to  the  target.  This  is  repeated  in  an 
iterative  process  using  updates  to  the  elasticity  parameters  of 
the  model  as  generated  by  a  multi-resolution,  non-linear 
optimization  scheme  in  order  to  achieve  a  suitable  match  in 
image  similarity.  Because  the  goal  of  the  reconstruction  is  to 


determine  a  spatial  mapping  of  tissue  elasticity,  this  process 
can  also  be  classified  as  an  inverse  problem. 

Our  observations  during  the  ongoing  development  and 
testing  of  this  method  have  prompted  questions  concerning 
the  quality  of  data  necessary  and  sufficient  to  achieve 
satisfactory  results  (i.e.  fidelity  of  the  reconstructed 
elasticity  image).  The  primary  inputs  to  the  reconstruction 
method  are  the  acquired  images  and  the  delineated  boundary 
conditions  on  the  region  of  interest.  While  it  is  clearly 
preferable  to  have  idealized  data,  in  reality,  both  inputs 
involve  varying  levels  of  manual  interaction.  As  an  initial 
study,  we  have  sought  to  test  the  effects  of  degradation  in 
data  quality  on  the  end  reconstruction  by  using  additive 
image  noise  and  randomized  boundary  condition  selection 
error. 

2.  METHODS 

2.1.  Elastographic  Reconstruction  Framework 

There  are  three  major  components  in  the  reconstruction 
framework:  a  biomechanical  model  of  tissue  response  to 
applied  deformation,  a  method  of  image  comparison,  and  an 
optimization  scheme.  For  the  current  version,  a  continuum- 
based  model  of  mechanical  equilibrium  using  isotropic 
Hookean  linear  elasticity  with  a  plane  stress  approximation 
is  employed  [5].  This  allows  for  a  reduction  of  the  general 
3D  formulation  of  Cauchy’s  Law  to  the  two  parameters  of 
Young’s  modulus  and  Poisson’s  ratio  in  2D.  The 
displacement  solution  of  the  finite  element  representation  of 
the  model,  solved  using  the  standard  Galerkin  method  of 
weighted  residuals  [6],  is  then  applied  to  the  nodes  of  a 
simple  triangular  mesh  based  on  the  source  image  domain  in 
order  to  perform  image  deformation.  The  mesh  is 
partitioned  by  K-means  clustering  (MATLAB  R14, 
Mathworks,  Nattuck,  MA)  into  N  number  of  regions,  each 
of  which  describes  a  distinct  set  of  homogeneous  elastic 
properties  for  a  grouping  of  adjacent  elements.  This  allows 
for  implementation  of  the  multi-resolution  approach  by 
creating  a  hierarchy  of  increasingly  finer  spatial 
distributions  of  elasticity  parameters,  which  has  been  shown 
to  be  an  improvement  upon  previous  versions  using  only  a 
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single  resolution  [2,3].  A  second  discretization  is  performed 
to  divide  the  target  image  into  M  number  of  rectangular 
zones  containing  approximately  equal  numbers  of  pixels. 
The  deformed  source  image  is  compared  to  the  target  using 
an  intensity-based  image  similarity  metric  (here,  the 
correlation  coefficient  [7])  in  the  evaluation  of  the  least 
squares  error  objective  function 

M 

(S  -S  V 

y^TRUE  ^EST) 

m=\ 

where  Strue  is  an  vector  of  the  (maximum)  similarity 
values  achieved  when  comparing  the  target  image  to  itself 
and  Sest  is  the  Mx7  vector  of  similarity  between  the  target 
and  model-deformed  source  image  created  using  current 
estimates  of  the  elastic  modulus  distribution.  It  should  be 
noted  that  Strue  has  by  definition  a  value  of  1  for  the 
correlation  coefficient. 

The  minimization  of  equation  (1)  using  a  Levenberg- 
Marquardt  approach  takes  the  form 

\j^J  +  al\/SE]=  -  Sest]  (2) 

where  J  is  the  Jacobian  matrix  of  size  MxN  estimating 
dS/dE,  AE  is  the  Nxl  vectors  of  updates  to  the  current 
elasticity  values,  and  a  is  the  scalar  regularization  term  for 
the  Hessian  matrix  as  described  in  [8]. 

2.2.  Material  Preparation  and  Image  Acquisition 

For  our  simulation  purposes,  a  two-material  skin  phantom 
had  been  previously  constructed  [2]  as  a  thin  membrane 
measuring  15  cm  x  15  cm,  with  a  single  5 -cm  circular  stiff 
inclusion  embedded  in  the  center  (Figure  1).  The  phantom 
was  manufactured  with  Smooth-On™  polyurethanes 
(Smooth-On,  Easton,  PA)  Evergreen  10  and  Evergreen  50. 
These  materials  have  essentially  indistinguishable  colors  but 
vary  significantly  in  their  elastic  modulus  values,  so  the 
former  was  used  as  the  bulk  material  and  the  latter  for  the 
inclusion.  Based  on  material  testing,  the  expected  contrast 
ratio  of  Young’s  modulus  values  was  determined  to  be 
approximately  5.7:1.  A  black  permanent  marker  was  used 
to  place  a  pattern  of  regularly  spaced  (-1  cm)  grid  lines  on 
the  membrane.  The  membrane  was  clamped  along  two 
opposite  edges  and  then  stretched  in  a  uniaxial  fashion  by 
means  of  a  milling  vise.  A  commercial  webcam  (Logitech 
QuickCam  Pro  4000)  was  mounted  above  the  assembly  to 
acquire  image  pairs  of  the  membrane  in  pre-  and  post- 
stretched  states  (960  x  1280  pixel  resolution,  8 -bit 
grayscale). 

2.3.  Reconstruction  Experiments 

Based  on  prior  work,  a  data  set  consisting  of  a 
particular  image  pair  and  associated  boundary  conditions 


Figure  1.  Experimental  phantom  membrane  system  (left)  and 
input  image  with  overlaid  finite  element  mesh  (right).  The 
inelusion  loeation  is  indieated  by  the  arrow  and  dotted  line.  The 
mesh  designates  the  aetual  region  reeonstrueted. 

known  to  produce  a  satisfactory  reconstruction  was 
designated  as  the  gold  standard  for  the  remainder  of  the 
experiments  (Figure  1).  In  order  to  test  the  effect  of 
increasing  amounts  of  additive  noise  on  the  reconstruction 
algorithm,  Gaussian  random  fields  of  1,  5,  10,  15,  20,  25, 
and  30%  noise  were  applied  to  the  base  target  image  in  three 
separate  trials.  This  presents  a  challenge  that  ascertains  the 
ability  of  the  similarity  metric  and  objective  function  to 
discern  a  proper  match. 

The  current  method  for  selecting  Dirichlet  boundary 
conditions  on  the  finite  element  mesh  is  semi-automated  and 
requires  the  user  to  make  a  final  determination  on  point 
correspondence.  The  second  experiment  was  intended  to 
simulate  the  targeting  error  of  the  user  (e.g.  visual  cues  and 
input  device  control).  Each  test  involved  applying 
randomized  vectors  of  equal  magnitude  to  alter  the 
boundary  conditions  of  the  gold  standard  data  set.  Errors  of 
0.1,  0.2,  0.3,  0.5,  0.75,  1.0,  1.5,  and  2.0  mesh  units  (scaled 
to  be  equivalent  to  pixel  coordinates)  were  used  in  two 
separate  trials  for  a  total  of  16  reconstructions.  Sub-pixel 
magnitudes  were  included  after  determining  that  the 
accuracy  of  selecting  a  feature  point  in  the  image/mesh  was 
typically  less  than  or  equal  to  0.5  units  for  users  ranging 
from  moderate  to  expert  skill. 

For  all  reconstructions,  resolutions  progressing  through 
TV  =  16,  36,  64,  144,  256,  and  400  regions  and  M  =  9 
similarity  zones  were  used;  domains  were  initialized  to 
homogeneous  elasticity  and  Poisson’s  ratio  held  constant  at 
0.485  to  represent  nearly  incompressible  mater ial(s). 

2.4.  Reconstruction  Analysis 

The  final  reconstructed  elasticity  values  were  modeled  as  a 
mixture  of  two  Gaussian  distributions,  and  a  threshold  was 
established  to  maximize  inter-class  variation  [9]  and 
subsequently  classify  each  region  as  bulk  or  stiff  material. 
Because  Dirichlet  boundary  conditions  are  exclusively  used 
in  these  reconstructions,  the  method  is  only  sensitive  to 
relative  differences  in  elasticity.  The  quantities  used  in 
evaluating  reconstruction  success  are  the  elasticity  contrast 
ratio,  localization  accuracy  of  the  inclusion,  and  an  overall 
measure  designated  the  ‘quality  of  reconstruction  score’ 
(QRS).  The  elasticity  contrast  ratio  (CR)  was  calculated 
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Figure  2.  Representative  reeonstmetions  with  image  noise.  From 
top  left:  1,  5,  10,  20,  25,  and  30%  additive  Gaussian  noise.  The 
reeonstruetions  are  visualized  as  two  materials,  with  whiter  areas 
indieating  higher  elastieity  eontrast  values. 


Table  1.  Reeonstruetion  quality  under  noise  eonditions 


Additive  image  noise 

%  Noise 

1 

5 

10 

15 

20 

25 

30 

LA 

0.92 

0.90 

0.91 

0.70 

0.69 

0.66 

0.56 

CR 

3.56 

3.45 

3.45 

3.24 

2.88 

2.83 

2.68 

Gold  standard:  LA  -  0.95,  CR  -  3.60 


Boundary  eondition  error 

Err 

0.1 

0.2 

0.3 

0.5 

0.75 

1.0 

1.5 

2.0 

AE 

0.96 

3.32 

2.21 

102 

0.93 

32.2 

12.6 

7.66 

LA 

0.87 

0.92 

0.88 

0.59 

0.94 

0.86 

0.86 

0.96 

CR 

3.63 

3.68 

3.44 

2.91 

3.46 

3.71 

3.78 

3.30 

CR  = 
AE  = 

elastieity  eontrast  ratio,  LA  =  loealization  aeeuraey 
initial  alignment  error  (%),  Err  =  error  magnitude. 

from  the  mean  values  of  the  two  material  classes,  and  the 
positive  predictive  value  of  identifying  stiff  material  within 
the  independently  segmented  boundary  of  the  inclusion 
gives  a  measure  of  localization  accuracy  (LA).  The  quality 
of  reconstruction  is  simply  then  the  product  QRS  =  CR*LA, 
which  allows  the  user  to  consider  the  other  two  measures  in 
conjunction. 


3.  RESULTS 

Figures  2  and  3  show  examples  of  reconstructions  achieved 
under  various  image  noise  and  boundary  condition  errors, 
and  individual  localization  errors  and  contrast  ratio  values 
are  listed  in  Table  1.  Note  that  the  data  for  the  image  noise 
experiment  was  averaged  from  the  three  trials,  and  that  the 
data  presented  for  the  boundary  condition  experiment  is 
from  one  [representative]  trial.  Figure  4  is  a  plot  of  the 
reconstruction  quality  decreasing  with  increasing  image 
noise,  and  Figure  5  shows  the  reconstruction  quality  trend 
plotted  against  the  change  in  initial  alignment  error  (detailed 
in  the  following  section)  relative  to  that  of  the  gold  standard. 


DDD 

ODD 


Figure  3.  Representative  reeonstmetions  with  boundary  eondition 
error.  Left  to  right:  0.1,  0.2,  0.3  units  (top  row);  0.75,  1.0,  2.0  units 
(middle  row,  trial  #1);  0.75,  1.0,  2.0  units  (bottom  row,  trial  #2). 
Error  magnitudes  greater  than  or  equal  to  0.5  mesh  units  are  not 
aeeurate  predietors  of  reeonstmetion  quality. 


4.  DISCUSSION 

From  visual  inspection  of  Figure  2,  it  is  apparent  that  the 
achieved  reconstruction  becomes  more  inaccurate  with 
increased  image  noise.  However,  the  ability  to  identify  and 
localize  the  stiff  inclusion  is  not  significantly  impaired  until 
a  noise  field  of  greater  than  10%  is  applied.  The  threshold 
was  found  by  determining  which  level  of  noise  provided  the 
best  minimum  sum  squared  error  fit  of  two  lines  to  the 
distribution  of  reconstruction  quality  in  Figure  4.  This 
would  indicate  that  the  similarity  metric  and  objective 
ftmction  are  robust  to  intensity  deviations  of  about  6  gray 
levels.  While  Gaussian  noise  is  one  of  several  possible  types 
and  may  not  always  be  an  ideal  model,  it  is  still  relevant  to 
acquisition  inaccuracy  and  corruption  processes  that  may  be 
encountered  across  several  medical  imaging  modalities.  The 
use  of  an  intensity-based  similarity  metric  appears  to  give 
the  method  an  advantage  in  being  generally  insensitive  to 
reasonably  expected  amounts  of  image  noise. 

Figure  3  demonstrates  that  because  of  the  random 
nature  of  the  boundary  condition  errors,  the  magnitude  is 
itself  not  an  accurate  indicator  of  reconstruction  quality. 
This  necessitated  the  introduction  of  a  more  suitable 
parameter  that  accounts  for  the  net  effect  of  the  altered 
boundary  conditions  in  order  to  perform  fair  evaluations.  In 
essence,  randomizing  the  vectors  at  every  node  causes  the 
optimization  to  use  an  unpredictable  starting  pose  and 
increases  its  chance  of  converging  to  an  improper  minimum. 
Therefore,  the  ‘initial  alignment  error’  (AE)  is  defined  as  the 
relative  percent  change  between  the  objective  function 
evaluation  using  the  gold  standard  boundary  conditions  and 
those  of  the  test  case.  An  as  example,  it  could  be  assumed 
that  vectors  of  magnitude  0.5  would  be  a  much  more 
tolerable  error  than  2.0,  but  it  is  the  significantly  larger  AE 
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Figure  4.  Reconstruction  quality  vs.  percent  additive  image  noise. 
The  drop-off  after  10%  additive  noise  indicates  the  threshold  of 
tolerance  for  the  method. 


Figure  5.  Reconstruction  quality  vs.  percent  change  in  initial 
alignment  relative  to  gold  standard.  The  majority  of  errors  tested 
produced  satisfactory  reconstructions. 


of  the  former  that  actually  predicts  the  poor  outcome. 
However,  it  should  also  be  noted  (results  not  shown  here) 
that  even  if  the  same  set  of  error  vectors  are  scaled  over 
varying  magnitudes,  there  is  no  clear  trend  in  alignment 
error  or  reconstruction  quality.  This  appears  to  imply  that 
certain  boundary  nodes,  most  likely  those  in  the  direction  of 
largest  strain,  have  a  greater  effect  on  reconstruction  and 
merit  particular  care  in  selection.  Other  factors  influencing 
unfavorable  reconstructions  are  most  likely  nonlinear  effects 
not  predicted  by  the  current  model  as  well  as  an  inherent 
lack  of  discrimination  by  intensity-based  similarity  metrics 
in  analyzing  the  regularity  of  the  imposed  grid  pattern. 
Nevertheless,  for  the  error  magnitudes  tested  that  best 
approximate  realistic  inaccuracies  (i.e.  <0.5  units),  the 
alignment  errors  were  small  and  quality  of  the  end 
reconstruction  was  seen  to  be  quite  good.  This  qualitatively 
validates  the  current  method  of  determining  point 
correspondence  as  a  reasonable  procedure  with  an 
accommodating  margin  (factor  of  four)  in  light  of  typical 
user  interaction. 

5.  CONCLUSIONS 

In  this  work,  we  have  presented  a  method  for  recovering 
elasticity  parameters  from  image  data  of  thin  membrane 
structures  by  maximizing  the  image  similarity  between  two 
different  states  of  mechanical  loading  within  the  context  of 
an  inverse  problem.  The  biomechanical  model,  multi¬ 
resolution  optimization,  and  image  acquisition  are  each 
modular  components  of  this  elastographic  reconstruction 
framework,  making  it  extensible  to  added  sophistication. 
Tests  were  conducted  to  examine  the  tolerance  of  the 
method  to  degraded  or  improper  inputs.  The  results  indicate 
that  the  gold  standard  data  set  was  mostly  optimal  for 
obtaining  a  successful  reconstruction.  Widening  disparities 
in  either  image  data  or  boundary  condition  selection  from 
that  in  the  gold  standard  caused  observable  trends  of 
declining  reconstruction  quality.  Based  on  these 
observations,  it  appears  that  the  method  handles  most 


expected  variations  encountered  in  image  acquisition  as  well 
as  the  majority  of  typical  user  inaccuracies.  Because  there 
are  complicated  effects  associated  with  mapping  of  the 
Dirichlet  boundary  conditions  in  constraining  the 
displacement  solution  of  the  model,  further  study  on  inter¬ 
rater  variability  in  selection  as  well  as  comparisons  with 
more  automated  point  correspondence  methods  is  ongoing. 
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BACKGROUND 

1 

RESULTS 

Changes  to  the  local  cytoarchitecture  induced  in  a  variety  of  pathoiogies  can  manifest 
as  aiterations  in  tissue  eiasticity  that  are  reievant  in  clinicai  examination  and  evaiuation. 
Many  eiastography  methods  are  typicaliy  dependent  on  the  specific  modaiity  around 
which  they  were  developed  (e.g.  magnetic  resonance  and  uitrasound  imaging).  We 
have  deveioped  ‘modality-independent  eiastography’  (MIE)  as  a  reconstruction  method 
that  recovers  the  materiai  properties  of  soft  tissue  via  model-based  analysis  of  image 
data  acquired  at  different  states  of  mechanical  loading.  The  algorithm  utilizes  image 
similarity  in  the  performance  of  a  multi-resolution,  non-linear  optimization.  Previous 
work  with  a  phantom  membrane  used  for  simulated  dermoscopic  applications  prompted 
this  preliminary  investigation  of  the  relative  effects  of  additive  image  noise  and 
boundary  condition  determination  errors  on  the  performance  of  the  method.  The  results 
as  quantified  by  elasticity  contrast  and  localization  accuracy  indicate  that  the 
reconstruction  process  is  robust  in  the  presence  of  realistic  levels  of  image  corruption 
and  tolerates  the  majority  of  boundary  condition  mapping  errors. 


PURPOSE 


The  inputs  to  the  reconstruction  process  are  in  two  major  forms:  image  data  and 
boundary  condition  estimation.  Inadequate  fidelity  in  either  quantity  is  capable  of 
affecting  the  success  of  the  reconstruction  through  some  form  of  model-data  mismatch. 
We  proposed  to  test  the  sensitivity  of  the  algorithm  to  various  levels  of  an  applied  noise 
process  by  altering  either  the  intensity  distribution  of  the  target  image  or  the 
displacement  vectors  defining  the  Dirichlet  boundary  conditions. 


METHODS 


Figure  1 .  Flow  chart  of  MIE.  After  acquisition,  source  and  target  images  (A)  are  discretized  into  regions  and  zones,  respectiveiy.  The  reconstruction  process 
invoives  updating  eiastic  moduius  vaiues  (B,E)  to  drive  a  finite  eiement  modei-based  image  deformation  (C)  untii  the  best  match  is  found  (D). 
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Figure  2.  MIE  reconstruction  experiment. 

(Left  panel)  A  two-material  phantom  mimicking 
skin  was  constructed  as  a  thin  membrane 
measuring  15  cm  x  15  cm,  with  a  single  5  cm 
circular  stiff  inclusion  embedded  in  the  center.  The 
phantom  was  manufactured  with  like-colored 
polyurethanes  which  have  an  inclusion-to-bulk 
elasticity  contrast  of  approximately  5.7:1.  The 
membrane  was  stretched  in  a  uniaxial  fashion 
while  a  CCD  camera  mounted  above  acquired 
image  pairs  of  the  membrane  in  pre-  and  post- 
deformed  states  (960  x  1280  pixel  resolution,  8-bit 
grayscale). 

(Right  panel)  Top  row:  Source  and  target  images 
with  overlay  of  finite  element  mesh  boundaries 
(red)  that  demarcate  the  area  reconstructed. 
Below:  Reconstruction  progression  over  increasing 
number  of  regions  (N  =  16,  64,  256,  400)  to  refine 
the  spatial  distribution  of  elasticity  values.  This 
reconstruction  serves  as  the  gold  standard  for  the 
remainder  of  this  work. 


Figure  3.  MIE  image  noise  reconstruction  experiment. 

Gaussian  random  fields  of  variable  strength  with  respect  to  the 
variance  of  non-background  pixel  values  were  applied  in  an 
additive  fashion  to  the  target  image.  Shown  in  the  left  column 
from  top  to  bottom  are  the  original  target  and  then  with  10%, 
20%,  and  30%  noise.  In  the  right  column  are  the  corresponding 
elasticity  reconstructions  after  application  of  a  thresholding 
scheme  to  classify  bulk  (black)  and  inclusion  materials 
(white/gray).  The  known  segmentation  of  the  inclusion  was  used 
to  retrospectively  calculate  the  positive  predictive  value  of 
identifying  the  correct  material  type  within  the  proper  boundaries 
as  well  as  the  mean  elasticity  contrast  of  the  overall  distribution. 
For  this  work,  our  overall  evaluation  of  reconstruction  quality  is 
expressed  as  the  product  of  these  two  quantities.  The  effect  of 
additive  noise  is  to  decrease  reconstruction  quality  as  evidenced 
in  the  progressively  poorer  localization  of  the  inclusion. 


Figure  4.  MIE  boundary  condition 


reconstruction  experiment. 
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Randomized  vectors  of  a  particular  magnitude  were  applied  to  the  boundary  condition  specifications  of  the  same  finite  element  mesh  used  for  all  reconstructions.  This 
simulates  targeting  error  by  the  user  in  the  currently  semi-automated  method  of  point  correspondence  selection,  and  the  effect  is  illustrated  in  the  top  row:  from  left  to 
right,  the  gold  standard  boundary  and  then  with  mis-estimation  of  0.75,  1.0,  and  2.0  mesh  units  (equivalent  to  pixel  coordinates)  in  the  Dirichlet  conditions  (slightly 
exaggerated  scale  for  visual  effect).  The  corresponding  reconstructions  in  the  middle  and  bottom  rows  demonstrate  that  two  different  trials  using  the  same  magnitude  of 
randomized  vectors  can  effect  very  different  levels  of  reconstruction  quality. 


Figure  5.  Reconstruction  quality  vs.  image  noise. 

Three  trials  of  image  noise  were  performed  (shown  averaged  and  with 
standard  error  bars);  the  drop-off  in  reconstruction  quality  indicates 
the  presence  of  a  threshold  at  approximately  1 0%  additive  Gaussian 

Figure  6.  Reconstruction  quaiity  vs.  boundary  condition 

Two  trials  of  eight  levels  of  noise  ranging  from  0.1  to  2.0  mesh  units 
were  performed.  Each  reconstruction  was  treated  as  a  separate  data 
point  based  on  its  initial  alignment  error,  defined  here  as  the  relative 
change  between  the  objective  function  evaluation  using  the  gold 
standard  boundary  conditions  and  those  of  the  [randomized]  test 


Figures  3  and  5  show  that  the  achieved  reconstruction  becomes  more  inaccurate  with  increased  image  noise.  However,  the  abiiity  to 
identify  and  locaiize  the  stiff  inclusion  is  not  significantiy  impaired  untii  a  noise  field  of  greater  than  10%  is  appiied.  The  threshoid  was 
found  by  determining  which  level  of  noise  provided  the  best  minimum  sum  squared  error  fit  of  two  iines  to  the  distribution  of 
reconstruction  quality.  This  wouid  indicate  that  the  simiiarity  metric  and  objective  function  are  robust  to  intensity  deviations  of  about  6 
gray  levels  in  an  8-bit  image.  While  Gaussian  noise  may  not  aiways  be  an  ideal  model,  as  a  preliminary  point  of  investigation,  it  is  stili 
relevant  to  acquisition  inaccuracy  and  corruption  processes  that  can  be  encountered  across  several  medical  imaging  modalities. 
Figure  4  demonstrates  that  the  magnitude  of  the  random  vectors  is  itseif  not  an  accurate  indicator  of  reconstruction  quaiity  because 
the  multipie  degrees  of  freedom  afforded  by  the  boundary  nodes  cause  the  optimization  to  use  an  unpredictabie  starting  pose, 
increasing  the  chances  of  converging  to  an  improper  iocai  minimum.  This  necessitated  the  introduction  of  the  initiai  aiignment  error 
(AE)  to  provide  a  consistent  means  of  comparison  between  triais  (Figure  6).  As  a  further  example,  it  couid  be  assumed  that  vectors 
of  magnitude  0.5  would  be  a  much  more  tolerable  error  than  2.0,  but  it  is  the  significantly  larger  AE  of  the  former  that  actually 
predicts  the  poor  outcome.  It  should  also  be  noted  (results  not  shown  here)  that  even  if  the  same  set  of  error  vectors  are  scaled  over 
varying  magnitudes,  there  is  no  clear  trend  in  alignment  error  or  reconstruction  quality.  This  appears  to  imply  that  certain  boundary 
nodes,  most  likely  those  in  the  direction  of  largest  strain,  have  a  greater  effect  on  reconstruction  and  merit  particular  care  in 
selection.  Nevertheless,  for  the  error  magnitudes  that  approximate  inaccuracy  in  boundary  condition  demarcation  (i.e.  <0.5  units),  the 
quality  of  those  reconstructions  was  satisfactory.  This  qualitatively  validates  the  current  method  of  determining  point  correspondence 
as  a  reasonable  procedure  with  an  accommodating  margin  (factor  of  four)  in  light  of  typical  user  interaction.  Further  research  is 
ongoing  into  validation  and  control  of  boundary  conditions,  as  well  as  more  automated  methods  of  point  correspondence. 
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The  use  of  palpation  information  for  skin  disease  characterization  is  not  as  commonly  used  as  in 
other  soft  tissues,  although  mechanical  differences  within  lesions  have  been  noted.  For  example, 
regions  of  hyperkeratosis  have  the  potential  to  transform  into  cancerous  lesions  and  likely  feature 
different  material  properties  from  those  of  surrounding  normal  tissue  due  to  varying  cytoarchitec- 
ture.  As  a  result,  the  spatial  distribution  of  lesion  mechanical  properties  may  serve  to  assist  a 
diagnosis  or  enhance  visualization  of  the  complete  extent  of  a  cancerous  region,  i.e.,  accurate 
information  regarding  the  margins  of  disease  for  surgical  therapy.  In  this  work,  a  multiresolution 
extension  to  a  novel  elastographic  imaging  method  called  Modality  Independent  Elastography 
(MIE)  is  used  to  characterize  the  mechanical  properties  of  a  skin-like  phantom  embedded  with  a 
mock  stiff  lesion.  Simulation  studies  were  also  performed  to  investigate  the  potential  for  charac¬ 
terizing  realistic  melanoma  lesions.  Elasticity  image  reconstructions  from  the  phantom  experiments 
localized  the  stiff  inclusion  and  had  good  correlation  between  the  Young’s  modulus  contrast  ratio 
and  experimental  measurements  from  material  testing.  In  addition,  multiresolution  MIE  was  shown 
to  be  a  more  robust  framework  than  its  single-resolution  version.  Results  from  the  melanoma 
simulation  demonstrate  the  potential  for  using  multiresolution  MIE  with  dermoscopic 
images.  ©  2005  American  Association  of  Physicists  in  Medicine.  [DOI:  10.1118/1.1895795] 


I.  INTRODUCTION 

Skin  cancers  are  a  growing  health  concern  in  the  United 
States,  with  total  annual  cases  being  reported  in  the  millions 
by  the  American  Cancer  Society.  There  are  three  major  types 
of  skin  cancers  [basal  cell  carcinoma  (BCC),  squamous  cell 
carcinoma  (SCC),  and  melanoma],  with  melanoma  estimated 
to  be  the  sixth  most  prevalent  cancer  and  an  estimated 
55,100  new  cases  (within  the  United  States)  to  be  diagnosed 
in  2004.^  In  general,  skin  cancers  develop  from  precancerous 
lesions  of  the  epidermis  that  have  dysplastic  changes  due  to 
the  damage  inflicted  by  ultraviolet  solar  radiation.  As  with 
other  cancers,  the  dysfunctional  cells  may  aggressively  com¬ 
pete  with  normal  tissue  for  nutrients  and  space.  The  progres¬ 
sion  from  a  benign  to  malignant  state  depends  upon  the  de¬ 
gree  of  cellular  differentiation  and  the  spatial  extent  of  the 
growth,  which  approximately  translates  into  the  pathological 
determination  of  grade  and  stage. 

When  skin  cancers  are  identihed  at  an  early  stage  and  are 
still  small  in  size,  surgical  excision  is  usually  straightforward 
and  effective.  If  the  disease  has  progressed  to  invade  deeper 
levels  of  the  skin,  treatment  becomes  more  difficult  and  may 
involve  more  invasive  surgery,  radiation,  and/or  chemo¬ 
therapy.  It  is  clear  that  the  early  detection  of  cancer  is  critical 
in  order  to  formulate  a  proper  treatment  plan  and  achieve  the 
most  favorable  clinical  outcome.  However,  detection  and  di¬ 
agnosis  still  rely  primarily  on  visual  inspection  followed  by  a 
biopsy  of  suspect  areas  for  histological  analysis.  Therefore,  a 
signihcant  proportion  of  diagnostic  technological  advances 
have  been  concerned  with  obtaining  a  better  view  of  the 
lesion  via  improved  optics  (i.e.,  dermoscopy)  or  more  ad¬ 
vanced  and  novel  imaging  systems  ranging  from  high- 
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frequency  ultrasound  to  confocal  laser  microscopy.  ’  Addi¬ 
tional  strategies  involving  electrical  impedance  mismatch,"^ 
Raman  spectroscopy,^  and  cytological  smears^  have  also 
been  forthcoming. 

As  opposed  to  other  methods  mentioned  above  which 
capitalize  on  electrical,  optical,  and  biochemical  phenomena, 
we  have  chosen  to  pursue  an  alternative  approach  to  skin 
health  assessment  which  is  based  on  its  mechanical  behavior. 
Detecting  changes  in  tissue  by  palpation  and  then  associating 
them  with  a  disease  state  has  had  a  longstanding  history  in 
clinical  medicine.  Although  a  health  assessment  of  skin  from 
palpation  is  performed  to  a  lesser  degree,  utilizing  changes  in 
the  mechanical  properties  to  characterize  the  skin  does  have 
precedent  within  clinical  dermatology.  One  thoughtful  re¬ 
view  by  Edwards  and  Marks  discusses  the  complex  mechani¬ 
cal  behavior  of  skin  when  subjected  to  in  vitro  and  in  vivo 
testing.^  Their  review  highlights  extensive  methodologies  be¬ 
ing  used  to  quantify  skin  mechanical  properties  (e.g., 
uniaxial  and  biaxial  extensometry,  torsion  stimulators,  inden- 
tometery,  bahistometric  tests,  shear  wave  application  de¬ 
vices,  dynamic  suction  methods,  ultrasonics,  and  electrody- 
namometry)  and  also  indicates  the  difficulties  in  comparing 
across  these  methods.  As  a  result,  Edwards  and  Marks  em¬ 
phasize  the  necessity  for  quantitative,  reproducible  methods 
to  assess  skin  health  given  the  wide  subjectivity  in  clinical 
analysis.^  Eor  example,  the  work  by  Draaijers  et  al  suggests 
that  reliable  subjective  assessment  of  the  pliability  of  scars 
requires  more  than  one  observer  while  measurements  using  a 
noninvasive  suction  device  can  be  accomplished  with  a 
single  observer.^  This  type  of  work  qualitatively  conhrms  the 
Edwards  and  Marks  conclusion  that  the  need  for  technology 
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and  automation  in  skin  assessment  will  be  essential  for  re¬ 
ducing  inter-rater  variability. 

While  the  characterization  of  skin  cancer  for  diagnostic 
purposes  and  possibly  surgical  intervention  is  an  interesting 
prospect,  other  investigations  have  begun  to  suggest  relation¬ 
ships  between  skin  elasticity  parameters  and  other  diseases. 
In  a  recent  study  using  a  noninvasive  suction  device,  Pierard 
et  al.  demonstrated  a  correlation  between  bone  mass  density 
(BMD)  and  skin  elasticity  parameters.  Specifically,  in  a  100- 
woman  study  in  which  a  portion  of  the  subjects  were  partici¬ 
pating  in  hormone  replacement  therapy,  a  positive  correlation 
existed  between  BMD  of  the  hip  and  femoral  neck  and  skin 
elasticity  parameters.  The  authors  clearly  state  that  their  goal 
was  not  to  develop  a  surrogate  BMD  assessment  test,  but  the 
results  are  nevertheless  intriguing. Using  a  similar  device, 
Yoon  et  al.  demonstrated  a  relationship  between  skin  elastic¬ 
ity  parameters  and  patients  afflicted  with  diabetes  mellitus.^^ 
Other  work  has  been  forthcoming that  demonstrates  the 
potential  for  using  noninvasive  measurements  of  skin  me¬ 
chanical  parameters  as  diagnostic  information. 

To  this  end,  the  field  of  elastography  has  established 
methods  to  spatially  characterize  the  mechanical  properties 
of  tissues  under  various  states  of  deformation  with  the  goal 
of  developing  functional  parameters  to  characterize 
disease.  ’  In  skin  cancer,  increases  in  cell  density,  atypia  in 
the  morphology  and  orientation  of  cells,  and  compositional 
alterations  (e.g.,  hyperkeratosis)  contribute  to  changes  in  the 
local  cytoarchitecture.  These  changes  in  mechanical  structure 
can  propagate  from  microscopic  to  macroscopic  levels  and 
may  manifest  as  a  distortion  of  the  normal  anatomy.  Given 
the  influence  of  mechanical  structure  on  the  behavior  of  de¬ 
forming  tissue,  elastographic  imaging  methods  may  be  well 
suited  for  detecting  and  monitoring  the  growth  of  these  can¬ 
cerous  anomalies.  In  fact,  advances  in  applying  ultrasound 
elastography  and  sonography  techniques  to  skin  are  being 
reported.  ’  Most  recently,  Gennisson  et  al  demonstrated 
the  use  of  a  new  sonoelastographic  probe  that  measured  a 
distinct  difference  between  dermis  and  hypodermis  shear 
wave  velocities  which  was  subsequently  used  to  estimate 
Young’s  modulus.  Although  interesting,  this  work  is  not 
completely  applicable  to  the  clinical  goals  of  understanding 
the  spatial  extents  of  a  melanoma  lesion. 

Following  previous  work  in  Ref.  23,  we  are  using  a  new 
elastographic  method  we  have  termed  “modality- 
independent  elastography”  (MIE)  that  combines  nonrigid  im¬ 
age  registration  with  an  elasticity  inverse  problem.  More  spe¬ 
cifically,  image  similarity  metrics  routinely  used  with  image 
registration  methods  are  recast  within  a  nonlinear  optimiza¬ 
tion  algorithm  whereby  mechanical  properties  (e.g..  Young’s 
modulus)  within  a  biomechanical  model  of  the  deforming 
tissue  become  the  driving  parameters  for  improved  image 
registration.  In  this  way,  the  MIE  method  circumvents  two 
potential  limitations  of  current  elastographic  techniques. 
Eirst,  it  is  not  inherently  dependent  on  preprocessing  steps 
such  as  homologous  feature  selection  and  tracking  which 
drive  active  contour  models  or  other  traditional 
displacement-based  iterative  methods^"^”^^  (however,  it  does 
require  the  determination  of  boundary  conditions).  Second, 
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because  it  is  an  image  processing  technique,  MIE  is  not  re¬ 
liant  on  a  particular  imaging  modality  such  as  in  ultrasound 
and  magnetic  resonance  elastography,  as  long  as  the  acquired 
images  provide  a  sufficient  pattern  to  allow  for  registration. 
Building  on  recently  completed  work  with  a  dual-mesh 
implementation,  in  this  paper  we  present  a  simplified  mul¬ 
tiresolution  elasticity  imaging  framework  for  Young’s  modu¬ 
lus  reconstruction.  In  addition,  phantom  and  simulation  ex¬ 
periments  demonstrate  its  utility  as  a  dermoscopic  image 
analysis  tool  for  evaluating  skin  lesions  based  on  material 
elasticity. 

As  a  final  point,  the  work  presented  here  represents  a 
potentially  new  application  of  the  MIE  approach  for  the 
characterization  of  skin  lesions  using  optical  images.  This 
may  have  significant  implications  at  many  length  scales 
(subcellular,  cellular,  matrix  level,  and  gross  tissue).  Eor  ex¬ 
ample,  properly  designed,  optically  based  MIE  could  be  used 
to  characterize  the  structural  development  of  tissues  at  the 
cellular  scale.  This  could  be  important  for  therapies  such  as 
Mohs  micrographic  surgery.  Mohs  is  a  surgical  technique 
which  combines  surgery  and  pathological  investigation  to 
more  effectively  remove  skin  tumors.  More  specifically,  after 
removing  visibly  cancerous  regions,  the  surgeon  removes  an 
additional  thin  layer  of  the  site  margin  and  creates  a  “map” 
of  the  border.  Upon  pathological  examination  of  the  removed 
layer,  the  “map”  can  be  used  to  target  the  remaining  cancer¬ 
ous  cells.  Currently,  the  Mohs  technique  is  a  time-consuming 
procedure,  but  the  success  of  the  procedure  is  compelling 
and  has  been  shown  to  be  cost  effective  with  certain 
considerations.  If  MIE  skin  imaging  could  accurately  assist 
or  replace  the  pathologic  characterization  of  the  margin  in 
less  time,  this  would  be  of  great  value  for  this  surgical 
therapy. 

II.  METHODS 

A.  Model  of  phantom/skin  elasticity 

One  critical  component  within  all  model-based  inverse 
problem  frameworks  is  the  selection  of  a  computational 
model  to  represent  the  continuum  of  interest.  In  our  phantom 
and  simulation  studies,  we  have  elected  to  employ  a  linear 
elastic  model  to  simulate  the  skin.  These  assumptions  (e.g., 
symmetry,  isotropy,  etc.)  allow  the  simplification  of 
Cauchy’s  law  from  36  stiffness  constants  to  2  and  employ  the 
equation 

V-o-=0,  (1) 

where  a  is  the  two-dimensional  (2-D)  Cartesian  stress  tensor 
and  is  defined  as 


T 
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The  constitutive  relationships  for  the  material  can  be  written 
as 
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where  E  is  the  Young’s  modulus,  v  is  Poisson’s  ratio,  and  u, 
V  are  displacements  in  the  x  and  directions,  respectively. 
For  this  work,  Poisson’s  ratio  was  assumed  to  be  constant  at 
0.485  for  our  skin  phantoms  and  tissue  simulations.  This 
value  was  found  by  searching  the  reconstruction  parameter 
space  for  an  optimal  value  that  achieved  maximum  similarity 
when  comparing  the  homogeneous  model-deformed  image 
to  its  acquired  counterpart.  The  constitutive  relationships  ex¬ 
pressed  in  (3)  represent  a  two-dimensional  approximation  to 
a  three-dimensional  system  which  assumes  a  symmetric,  iso¬ 
tropic,  thin  specimen  in  equilibrium  and  stresses  that  are 
constrained  to  lie  within  the  plane,  i.e.  the  classic  plane  stress 
approximation.^^  Using  the  Galerkin  method  of  weighted  re¬ 
siduals  to  integrate  this  set  of  partial  differential  equations,  a 
finite  element  framework  is  generated  and  can  be  solved  to 
represent  a  displacement  field  for  a  given  distribution  of 
Young’s  modulus  values. The  boundary  conditions  for  our 
studies  below  were  either  manually  derived  from  a  structured 
grid  representation  as  in  the  phantom  system  or  prescribed 
by  the  user  in  the  case  of  the  simulation  studies. 

B.  Modality-independent  eiastography  (MIE) 

The  MIE  framework  begins  with  the  acquisition  of  a 
baseline  predeformed  “source”  image  and  a  post-deformed 
“target”  image.  The  “source”  image  set  is  used  to  create  a 
well-resolved  finite  element  mesh  of  the  tissue  domain.  In 
previous  work,  a  second  coarse  mesh  was  also  specified  on 
the  domain  and  was  used  specifically  as  the  mechanical 
property  reconstruction  grid.  In  this  work,  a  new  single¬ 
mesh  region-based  multiresolution  MIE  approach  has  been 
employed  which  simplifies  previous  dual-grid  techniques 
with  the  generation  of  a  structured  regionalization  using  a 
^-means  clustering  algorithm  based  on  the  element  centroids 
of  the  well-resolved  mesh.  A  ^-means  clustering  algorithm 
iteratively  partitions  the  element  centroids  into  a  given  num¬ 
ber  (K)  of  regions  (where  K  is  the  user-defined  number  of 
desired  clusters)  such  that  the  sum  of  all  point-to-r^g/6>^  cen¬ 
troid  distances  over  all  regions  is  minimized.  The  advantage 
of  using  the  ^-means  clustering  approach  as  opposed  to  a 
regular  grid  is  that  the  clustering  approach  can  more  appro¬ 
priately  fit  irregular  domains  (e.g.  the  circular  domain  for  the 
dermoscopic  image  set).  Eor  this  work,  the  implementation 
in  the  MATLAB  (MathWorks,  Natick,  MA — 
www.mathworks.com)  statistics  toolbox  was  used.  Eigure  1 
illustrates  an  example  of  this  approach  on  a  circular  domain 
whereby  the  element  centroids  have  been  clustered  into  16 
separate  homogeneous,  isotropic  regions. 

The  method  has  been  adapted  to  a  multiresolution  strategy 
whereby  coarser  resolutions  (i.e.,  fewer  regions)  can  be  ini- 
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Fig.  1.  .^T-means  material  property  clustering  for  a  circular  domain  with  16 
property  regions  designated. 


tially  reconstructed  to  provide  better  initial  guesses  to  subse¬ 
quent  resolutions.  The  use  of  hierarchical  multiresolution 
structures  within  both  rigid  and  nonrigid  registration  algo¬ 
rithms  has  a  longstanding  precedent  and  lends  credence  to  its 
application  here.^^“^^  In  this  work,  six  progressively  finer 
resolutions  were  used  within  each  reconstruction  (16,  36,  64, 
144,  256,  400  regions). 

Once  the  mesh  and  ^-means  resolutions  have  been  speci¬ 
fied,  the  reconstruction  algorithm  begins  by  assigning  an  ini¬ 
tial  modulus  value  to  each  region  (a  homogeneous  initializa¬ 
tion  is  assumed)  at  the  first  resolution,  weighted  residual 
equations  are  integrated,  boundary  conditions  are  applied, 
and  the  matrix  equation  system  is  generated: 

VA{EMu}  =  {b},  (4) 

where  [A(£’^)]  represents  the  model  stiffness  matrix  based  on 
the  current  distribution  of  properties  {u}  is  the  vector  of 
unknown  tissue  displacements,  and  {b}  is  the  vector  of 
known  forces  acting  on  the  system  and  boundary  conditions. 
Upon  the  calculation  of  tissue  displacements,  the  source  im¬ 
age  can  be  deformed.  This  model-deformed  source  image  is 
then  compared  to  the  target  image  using  an  image  similarity 
method  ’  which  is  calculated  over  a  number  of  discrete 
spatial  zones  (e.g.,  for  all  reconstructions,  approximately  400 
similarity  zones  were  designated  within  the  image  for  a  com¬ 
parison).  Modulus  values  in  the  regions  are  updated  based  on 
maximizing  the  similarity  between  the  deformed  source  im¬ 
age  and  the  target  image  over  all  the  similarity  zones  until  a 
tolerance  is  reached  or  the  desired  number  of  iterations  has 
been  completed. 

With  respect  to  the  optimization  framework  for  MIE,  it 
can  be  represented  as  a  least  squared  error  objective  func¬ 
tion: 

</.(£)  =  min{||5(i,)-5(4)|P},  (5) 

where  S{E^  is  the  similarity  value  achieved  when  comparing 
the  target  image  to  itself  (i.e.,  the  maximum  value  for  the 
similarity  metric)  and  S{E^  is  the  similarity  between  the 
model-deformed  source  image  and  the  target  image  using  the 
current  estimate  of  the  elastic  modulus,  E^.  Equation  (5)  can 
be  solved  by  employing  a  Newton-Raphson-based  approach: 
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Fig.  2.  Multiresolution  MIE  algorithm 
flow  chart  where  “/?  =  1 , 2 , 3  •  •  • ,  ^max ’’ 
is  the  resolution  level  with  Rmax  l^e 
most  well  resolved;  and  is  the 
number  of  material  regions  within  a 
particular  resolution 


[[7^[7]  +  a[/]]{Ai}  =  -  A(4)},  (6) 

where  [/]  is  the  MXN  Jacobian  matrix  of  the  form  J 
=  dS{E^I  dE,  M  is  the  number  of  similarity  measurement 
zones,  and  N  is  the  number  of  material  property  regions  and 
is  equivalent  to  K  as  designated  in  the  ^-means  clustering 
algorithm.  The  details  of  Eq.  (6)  have  been  reported 
previously.^^’^^  Because  [/^][/]  (an  approximation  to  the 
Hessian  matrix)  tends  to  be  ill  conditioned,  it  is  regularized 


with  an  empirically  determined  a  parameter  found  in  the 
standard  Levenberg-Marquardt  approach.^^  The  determina¬ 
tion  of  this  regularization  parameter  is  described  in  Ref.  37. 
Figure  2  is  a  flow  chart  of  the  new  multiresolution  MIE 
approach. 

In  previous  work,  we  have  analyzed  the  performance  of 
our  MIE  algorithm  with  respect  to  four  standard  image  simi¬ 
larity  metrics  found  within  the  literature:  the  sum  of  squared 
differences,  normalized  mutual  information,  the  correlation 
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Fig.  3.  Experimental  data  from  the  skin- stretching  setup  shown  in  Fig.  4:  (a)  baseline,  (b)  5  mm,  (c)  10  mm,  (d)  15  mm,  (e)  20  mm. 


coefficient  (CC),  and  the  gradient  correlation  coefficient 
(GC).^^  Within  this  work,  the  correlation  coefficient  and  gra¬ 
dient  correlation  coefficient  were  used  for  the  similarity  mea¬ 
surements. 

Briefly  stated,  the  CC  operates  on  the  distribution  and 
mean  intensity  values  of  the  overlapping  regions  of  two  im¬ 
ages  where  /;  would  represent  the  intensity  values  within  the 
“target”  image  and  I2  would  be  the  model-deformed  “source” 
image.  The  correlation  coefficient  can  be  calculated  by  the 
expression 


s,ai(0-/i)(/20')-/2) 

hi(h{l)-h)\l2{l)-l2? 


V  G  /i  n  4, 


(7) 


where  /j,  I2  are  the  mean  intensity  values  within  each  respec¬ 
tive  image,  and  i  is  the  ith  pixel  within  the  respective  image. 
The  CC  metric  is  calculated  by  applying  the  correlation  co¬ 
efficient  to  images  that  have  been  processed  by  any  of  the 
standard  edge  detection  functions  (e.g..  Canny,  Sobel,  etc.). 


C.  Phantom  construction 

A  phantom  was  constructed  that  was  approximately  25  cm 
long,  15  cm  wide,  and  approximately  2  mm  thick.  The 
inclusion-surrounding  bulk  material  of  the  phantom  was 
Smooth-On^M  Evergreen  10  polyurethane  with  an  additive  to 
allow  permanent  marker  to  adhere  to  the  material  surface 
(Smooth-On,  2000  Saint  John  Street,  Easton,  PA).  A  cylin¬ 
drical  inclusion  was  placed  centrally  within  the  membrane 
phantom  that  was  approximately  5  cm  in  diameter  and  was 
made  of  a  stiffer  polyurethane  material  (Smooth-On^M  Ever¬ 
green  50).  The  inclusion  material  was  chosen  for  its  relative 
stiffness  to  that  of  Evergreen  10  and  its  color  which  is  the 
same  (to  study  the  case  of  non-pigmented  lesions).  After  the 
phantoms  had  set,  a  permanent  marker  was  used  to  draw 


15  cmX  15  cm  grid  with  1  cmX  1  cm  squares  on  the  phan¬ 
tom  surface.  Eigure  3(a)  shows  the  skin  phantom  used  for 
data  collection  in  this  series  of  experiments. 


D.  Image  acquisition  protocoi 

To  acquire  the  pre-  and  post-deformed  images  of  the 
stretched  skin  phantom,  the  membrane  was  first  secured  in 
customized  clamps  attached  to  a  milling  vise  to  form  a  trans¬ 
lation  stage  and  then  brought  level  with  a  nominal  applied 
load  to  define  the  baseline  position.  Images  were  taken  by  a 
commercial  web  camera  (Logitech  QuickCam  Pro  4000, 
960  X  1280  pixel  resolution)  that  was  rigidly  mounted  above 
the  membrane  at  a  single  location  to  ensure  a  fixed  field  of 
view  and  frame  of  reference  for  the  duration  of  the  experi¬ 
ment.  A  series  of  five  total  images  was  collected  (eight-bit 
grayscale)  via  laptop  control  of  the  camera — the  baseline 
predeformation  position  and  four  subsequent  positions  with 
incremental  stretches  of  approximately  5  mm  each.  Eigure  4 


Fig.  4.  An  illustration  of  the  skin-phantom  setup  for  image  acquisition. 
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Fig.  5.  (a)  Melanoma  lesion,  reproduced  with  the  permission  of  Dr.  Lehmann,  M.D.,  ©  Dermatlas,  www.dermatlas.org,  melanoma_l_040510.  (b)  Simulated 
horizontal  and  (c)  vertical  displacements  shown  (axis  references  are  in  meters  while  the  gray  scale  is  in  microns).  It  should  be  noted  that  the  contained  region 
within  the  border  represents  the  spatial  regions  of  stiffness  in  this  simulation  and  was  not  contained  within  image  data  provided  to  the  MIE  algorithm. 


is  a  schematic  of  the  experimental  setup,  while  Fig.  3(a)-3(e) 
shows  an  example  dataset. 

E.  Material  testing  protocol 

Material  testing  was  performed  in  order  to  determine  the 
accuracy  of  the  reconstructed  Young’s  modulus  values. 
When  the  phantoms  were  poured,  specimens  of  both  the  bulk 
and  stiff  polyurethane  were  allowed  to  cure  in  separate  con¬ 
tainers  from  the  membranes.  These  samples  were  then  cut 
into  1  cm  XI  cm  XI  cm  cubes.  Compression  testing  was 
performed  on  an  EnduraTEC  ELE  3200  material  tester  (En- 
duraTEC  Systems  Group,  Minnetonka,  MN).  The  polyure¬ 
thane  was  assumed  to  be  elastic,  homogenous,  and  isotropic. 

The  Enduratec  material  testing  protocol  involved  ramping 
the  actuator  linearly  from  the  zero  position  to  24%  strain  at 
2%  strain  increments.  The  max  strain  value  was  chosen  to 
extend  slightly  beyond  the  range  of  observed  strain  in  the 
experiment  shown  in  Eig.  3  which  was  approximately  22% 
strain  for  the  bulk  material.  Although  the  stiffer  inclusion 
material  only  experienced  approximately  l%-2%  strain, 
stress  values  were  reported  for  the  full  strain  range  up  to 
24%.  Between  each  change  in  axial  position  a  three  second 
dwell  time  was  imposed  to  allow  viscoelastic  responses  to 
subside.  Stress-strain  plots  were  produced  for  both  the  bulk 
material  and  the  inclusion  material.  Three  samples  of  each 
material  were  tested  and  an  average  curve  was  calculated. 

R  Phantom  experiment 

To  quantitatively  test  the  MIE  method  within  the  context 
of  dermoscopic  applications  using  optical  images,  a  series  of 
studies  using  the  elastic  membrane  of  Eig.  3  were  employed 
within  the  setup  of  Eig.  4.  The  single  inclusion  phantom  was 
considered  to  be  representative  of  a  single  lesion  on  the  skin 
surface  (nonpigmented  in  this  case).  The  multiresolution 
MIE  technique  was  used  at  each  successive  deformation  for 
a  total  of  four  elasticity  image  reconstructions  per  similarity 
metric  (in  this  case  CC  and  GC  only).  The  computational 
domain  for  these  calculations  involved  1051  nodes  and  1973 
elements.  Boundary  conditions  were  generated  by  analyzing 


the  pre-  and  post-deformed  structured  grid  and  estimating  the 
domain’s  deformation.  The  Young’s  modulus  reconstructions 
were  then  compared  to  the  elasticity  values  as  generated 
from  the  material  testing  protocol.  It  should  be  noted  that 
only  Young’s  modulus  contrast  was  compared  in  these  evalu¬ 
ations.  This  is  due  to  the  manner  in  which  boundary  condi¬ 
tions  are  prescribed  in  the  model  system.  Currently,  the  ap¬ 
proach  is  driven  by  displacement  boundary  conditions  (i.e., 
Dirichlet  type)  which  consequently  make  the  elastic  model 
only  sensitive  to  Young’s  modulus  contrast.  Without  knowl¬ 
edge  of  an  applied  stress  at  the  boundary  or  a  prescribed 
material  property  within  the  domain,  absolute  properties  can¬ 
not  be  determined.  In  addition,  it  must  also  be  noted  that  the 
reconstructions  were  constrained  to  a  region  of  the  phantom 
that  was  smaller  than  the  overall  phantom.  This  was  a  result 
from  observing  that  at  higher  stretch  states,  out-of-plane  dis¬ 
tortions  of  the  membrane  became  more  prominent  in  the  pe¬ 
riphery. 

G.  Simulation  experiments 

In  an  effort  to  test  the  algorithm  within  the  context  of  a 
more  realistic  image  acquisition  scenario  for  skin  cancer,  a 
simulation  study  was  performed  on  an  image  of  the  1  cm 
melanoma  lesion  shown  in  Eig.  5(a).  In  addition,  a  grid  struc¬ 
ture  was  not  specifically  applied  to  the  lesion  image  so  as  to 
test  whether  the  natural  skin-texture  itself  contained  suffi¬ 
cient  image  information  for  reconstruction.  The  lesion  was 
provided  by  the  Dermatology  Image  Atlas  project  (www.der- 
matlas.org.  Image  Name:  melanoma_I_040510,  Contributed 
by  Eric  Ehrsam,  M.D.)  and  represents  a  1  cm  pigmented 
melanoma  plaque,  located  on  the  left  arm  of  a  35-year-old 
woman.  Eor  the  simulation  boundary  conditions,  an 
annulus-shaped  mechanical  stretching  device  was  assumed 
which  would  systematically  stretch  two  semicircular  regions 
apart  by  2  mm.  The  melanoma  was  assumed  to  have  a  2:1 
elasticity  contrast  level  with  normal  tissue,  i.e.,  the  mela¬ 
noma  was  twice  as  stiff  as  the  surrounding  skin.  The  com¬ 
putational  domain  for  the  inverse  problem  contained  1294 
nodes  and  2459  elements.  In  addition,  the  mesh  used  to  gen¬ 
erate  the  forward-problem  data  was  approximately  10% 
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various  levels  of  bulk  material  strain  (for  all  ratio  determinations,  the  inclusion’s  2%  strain  value  for  Young’s  modulus  was  used  which  was  the  approximate 
maximum  strain  reached  in  the  inclusion  based  on  experiment  observations). 


more  resolved  than  the  mesh  used  for  reconstruction.  This 
introduced  a  small  degree  of  variability  to  the  boundary  con¬ 
ditions  and  image  deformation  to  simulate  potential  acquisi¬ 
tion  noise.  Figure  5(b)  and  Fig.  5(c)  illustrates  the  localized 
application  of  the  stretching  and  the  simulated  solution  for 
both  horizontal  and  vertical  displacements,  respectively. 
Upon  completion,  these  image  data  was  used  as  input  to  the 
multiresolution  MIE  algorithm.  Results  are  reported  using 
the  CC  and  GC  image  similarity  methods. 

III.  RESULTS 
A.  Material  testing 

During  the  material  testing  phase,  additional  cyclic  testing 
was  performed  in  which  viscoelastic  behavior  was  noted.  As 
a  result,  a  waiting  period  was  utilized  at  each  strain  level  to 
allow  viscoelastic  responses  to  subside.  The  stress/strain  be¬ 
haviors  at  these  quasistatic  time  periods  for  the  bulk  material 
and  inclusion  are  shown  in  Figs.  6(a)  and  6(b).  In  addition. 


discrete  finite  difference  approximations  were  made  at  the 
various  strain  levels  to  estimate  a  Young’s  modulus  value. 
Table  I  represents  a  distribution  of  values  calculated  within 
the  strain  ranges  tested.  Once  each  modulus  was  calculated 
for  each  acquired  strain  level,  a  distribution  of  Young’s 
modulus  contrast  ratios  was  calculated  and  is  expressed  by 
Eq.  (2): 


Table  I.  Young’s  elastic  modulus  values  at  several  experimental  strain  lev¬ 
els. 


Material  strain  (%) 

Bulk  material  (kPa) 

Inclusion  material  (kPa) 

2 

112.5 

868.0 

8 

154.3 

1635.0 

12 

180.0 

2125.0 

16 

182.5 

2375.0 

23 

206.3 

2412.5 
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Table  IL  Contrast  ratios  for  each  experimental  strain  level. 


Bulk  material  strain 

Contrast  ratio 

Stretch  1:  8% 

5.7 

Stretch  2:  12% 

5.0 

Stretch  3:  16% 

4.6 

Stretch  4:  23% 

4.2 

CR(sbuik)=  (8) 

^bulk 

whereby  the  inclusion’s  2%  strain  value  (approximate  maxi¬ 
mum  strain  reached  in  the  inclusion  based  on  experimental 
observations)  for  Young’s  modulus  was  used  and  the  bulk 
material  was  allowed  to  vary  over  the  entire  strain  range. 
This  contrast  ratio  formulation  reflects  the  reality  of  the 
membrane  experiments  shown  in  Fig.  3  whereby  the  soft 
surrounding  material  experienced  the  majority  of  deforma¬ 
tion  with  the  inclusion  remaining  relatively  unchanged  over 
all  applications  of  stretch.  The  distribution  of  the  contrast 
ratios  as  described  by  Eq.  (8)  at  differing  strain  levels  is 
shown  in  Fig.  6(c). 


31 

To  assist  in  determining  inclusion-to-bulk  contrast  ratios 
for  different  bulk  strain  levels  in  each  experimental  stretch  as 
reflected  in  Fig.  3,  an  exponential  curve  fit  was  prescribed: 

CR(8  =  Sb„,k)=A  +  B^>-c^  (9) 

whereby  A =4.0,  .S=5.0,  C=13.8.  The  root  mean  square  con¬ 
trast  ratio  error  between  model  and  acquired  data  over  the 
entire  acquired  strain  range  was  0.093.  The  quality  of  the 
exponential  model  can  be  seen  in  Fig.  6(c).  Using  the  expres¬ 
sion  described  in  (9),  a  series  of  Young’s  modulus  contrast 
ratios  values  could  be  tabulated  as  a  function  of  the  specific 
strain  levels  used  within  the  experiments  of  Fig.  1.  These 
levels  were  determined  by  manually  measuring  strain  levels 
within  regions  of  the  bulk  material  from  the  optical  images. 
Table  II  reports  the  approximate  contrast  ratio  for  Young’s 
modulus  at  the  various  bulk  material  strain  levels  experience 
during  the  stretching  experiments  using  Eq.  (9). 

B.  Multiresolution  MIE  phantom  reconstructions 

Figures  7  and  8  are  representations  of  the  multiresolution 
elasticity  image  reconstruction  performance  for  each  of  the 
different  stretch  states  shown  in  Fig.  3  using  CC  and  GC  as 
the  basis  for  image  similarity,  respectively.  The  boundary 


Normalized  Transect  Distance 


Normalized  Transect  Distance  Normalized  Transect  Distance 


Normalized  Transect  Distance 


3a-3b  3a-3c  3a-3d  3a-3e 


Fig.  7.  An  illustration  of  elasticity  image  reconstructions  using  CC  where  each  column  represents  the  respective  stretch  relative  to  Fig.  3  (e.g.,  3a-3b 
represents  the  stretch  from  base  to  the  first  increment).  The  top  image  shows  the  location  of  a  transect  as  designated  by  the  T  and  the  gray  scale  associated 
with  Young’s  Modulus  (Pa).  The  middle  row  represents  the  reconstructed  elasticity  images  at  each  stretch  state  with  the  contour  showing  the  actual  inclusion 
location.  The  bottom  row  shows  the  elastic  property  contrast  ratio  as  compared  to  that  predicted  with  the  material  testing  (shown  as  dark  box-like  contour) 
along  the  transect  T. 
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Fig.  8.  Illustration  of  elasticity  image  reconstructions  using  GC  where  each  column  represents  the  respective  stretch  relative  to  Fig.  3  (e.g.,  3a-3b  represents 
the  stretch  from  the  base  to  the  first  increment).  The  top  row  represents  the  reconstructed  elasticity  images  at  each  stretch  state.  The  bottom  row  shows  the 
elastic  property  contrast  ratio  as  compared  to  that  predicted  with  the  material  testing  (shown  as  a  dark  box-like  contour)  along  the  transect  T  which  was 
designated  in  Fig.  7. 


contour  represents  the  inclusion  location  as  shown  within  the 
image.  The  contrast  ratios  designated  within  the  transect  im¬ 
ages  of  Figs.  7  and  8  were  based  on  Table  II. 

Figure  9  shows  elasticity  images  at  varying  stage  resolu¬ 
tions  within  the  multiresolution  MIE  reconstruction  (recon¬ 
struction  shown  is  the  GC-3a-3b  stretch  regime). 

To  test  the  effect  of  the  multiresolution  framework,  the 
same  optical  images  were  provided  to  our  algorithm  using  a 
400  property  region  resolution  with  an  initial  guess  of  homo¬ 
geneity  (i.e.,  coarser  resolution  solutions  were  not  used  as 
initial  guesses).  In  results  not  reported  here,  the  CC  recon¬ 
struction  was  able  to  localize  and  quantify  the  stiff  region 
similar  to  that  of  Fig.  7  at  the  high  stretch  states  but  was 
much  worse  with  respect  to  the  initial  stretch  state  (i.e., 
3a-3b  image  reconstruction). 

Figure  10  represents  the  GC  result  for  the  four  stretch 
states  using  the  single  400  region  high  resolution  parametri- 
zation.  In  Fig.  10(a)  (3a-3b  stretch  state),  the  inclusion  is 


16  Regions  36  Regions  64  Regions 


144  Regions  256  Regions  400  Regions 


Fig.  9.  Elasticity  image  reconstruction  for  the  GC  3a-3b  reconstruction  case 
at  various  resolutions  of  the  multiresolution  algorithm.  The  gray  scale  is  the 
same  as  in  Fig.  7. 


localized  but  the  contrast  resolution  is  poor  compared  with 
its  multiresolution  counterpart  in  Fig.  8,  first  column.  At  sub¬ 
sequent  stretch  states  (3a-3c,  3a-3d,  and  3a-3e),  the  elastic¬ 
ity  image  has  not  converged  to  an  acceptable  representation 
of  the  inclusion.  Interestingly,  the  distance  traveled  by  grid 
squares  within  the  homogenous  regions  near  the  stretching 
edge  within  the  second  stretch  state  (3a-3c)  is  approximately 
the  size  of  one  grid  square  (~1  cm). 

It  is  evident  that  by  using  a  single  high  resolution  param- 
etrization  as  opposed  to  a  multiresolution  approach,  a  local 
minimum  is  found  and  the  elasticity  image  degrades  consid¬ 
erably.  Consequently,  the  error  magnitude  for  the  image 
shown  in  Fig.  8,  the  second  column  is  a  factor  of  50% 
smaller  than  that  of  Fig.  10(b)  thus  demonstrating  that  Fig. 
10(b)  indeed  represents  a  local  minimum  (it  should  be  noted 
that  all  parameters  were  identical — number  of  similarity 
zones,  filtering,  regularization,  relaxation,  etc.). 

C.  Multiresolution  MIE  melanoma  reconstruction 
simulations 

In  addition  to  the  experimental  results  shown  above,  sev¬ 
eral  similar  simulations  were  executed  using  a  pigmented 
melanoma  image.  Figure  11  shows  the  elasticity  image  re¬ 
construction  and  transect  results  using  the  multiresolution 
MIE  framework  for  both  CC  and  GC.  Eigure  12  illustrates 
the  inter-resolution  results  from  the  GC  reconstruction 
shown  in  Eig.  11. 

IV.  DISCUSSION 

The  elasticity  image  results  from  phantom  (Eigs.  7-9)  and 
simulation  (Eigs.  11  and  12)  studies  demonstrate  the  utility 
of  the  multiresolution  MIE  approach.  In  addition,  comparing 
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Fig.  10.  GC  reconstructions  using  single  400  property 
zone  resolution  for  (a)  3a-3b,  (b)  3a-3c,  (c)  3a-3d,  and 
(d)  3a-3e,  respectively.  The  gray  scale  is  the  same  as  in 
Fig.  7. 


the  results  in  Figs.  8  and  10  clearly  illustrates  that  instances 
can  exist  in  which  a  single-resolution  approach  will  fail 
whereas  the  mutliresolution  succeeds.  A  separate  but  related 
concern  which  is  still  under  investigation  is  the  degree  and 
content  of  the  image  pattern  needed  to  facilitate  reconstruc¬ 
tion;  however,  the  preliminary  elasticity  image  results  from 
the  melanoma  simulations  reported  herein  suggest  that  a  suf¬ 
ficient  intensity  content  exists  in  standard  dermoscopic  im¬ 
ages. 

Another  important  advance  in  this  paper  over  previous 
work  is  the  comparison  between  reconstructed  elastic  prop¬ 
erties  and  their  separately  measured  counterparts.  The  stress- 
strain  curves  shown  in  Figs.  6(a)  and  6(b)  and  modulus  val¬ 
ues  in  Table  I  demonstrate  a  nonlinear  elastic  behavior.  A 
good  representative  exponential  fit  to  the  Young’s  modulus 
contrast  ratio  data  was  achieved  in  Fig.  6(c)  and  provides  a 
direct  comparison  to  MIE-derived  Young’s  modulus  proper¬ 
ties.  One  shortcoming  is  that  because  MIE  is  completely 
driven  by  displacement  boundary  conditions,  only  the  con¬ 
trast  in  Young’s  modulus  values  can  be  compared.  However, 
the  goal  within  this  work  is  to  investigate  elastic  properties 
as  a  mechanism  for  contrast  within  medical  images. 

Overall,  the  elastic  image  reconstructions  shown  in  Eigs. 
7  and  8  demonstrated  good  localization  with  a  varied  perfor¬ 
mance  in  maintaining  lesion  shape  integrity  for  both  the  CC 
and  GC  similarity  methods,  respectively.  It  appears  that  at 
high  strain  levels,  MIE  was  less  successful  at  capturing  the 
anticipated  contrast  ratio.  In  fact,  in  both  CC  and  GC,  the 
ratio  was  overestimated,  thus  producing  more  contrast.  It 
should  be  noted  that  the  reconstructions  shown  were  per¬ 
formed  on  a  domain  that  represented  only  a  portion  of  the 
image  that  surrounded  the  inclusions  (—3-4  cm  from  the 
inclusion  border).  This  was  due  to  our  inability  to  completely 
control  the  physical  boundaries  of  the  phantom  given  the 
large  mismatch  between  the  stiffness  values  of  the  two  ma¬ 
terials.  This  manifested  itself  as  out-of-plane  warping  of  the 


phantoms,  i.e.,  a  wrinkling  at  edges  as  the  strain  on  the  skin 
phantoms  increased.  The  spatial  location  of  these  membrane 
distortions  was  more  prominent  with  the  distance  from  the 
inclusion.  By  making  a  more  localized  reconstruction  region, 
the  influence  of  these  distortions  was  minimized  although 
some  effects  are  undoubtedly  present.  Ultimately,  these  out- 
of-plane  motions  would  be  interpreted  as  planar  strains  in  the 
optical  image  acquisition  system  shown  in  Eig.  4.  Although 
this  variability  in  shape  integrity  existed,  successful  localiza¬ 
tion  was  achieved  for  all  stretch  states.  It  was  encouraging 
that  at  small  stretch  states,  where  the  model  is  most  appro¬ 
priate,  proper  quantitative  contrast  ratios  were  achieved 
(stretch  states  3a-b,  3a-c  in  Eigs.  7  and  8).  Eurther  encour¬ 
agement  was  provided  by  successful  localizations  at  high 
stretch  states  whereby  nonlinear  behavior  is  undoubtedly 
present  and  the  small-strain  assumptions  are  compromised 
(although  the  quantitative  contrast  ratio  was  not  as  satisfy¬ 
ing).  Undoubtedly,  a  large-deformation  model  is  necessary  at 
these  higher  strains  to  match  contrast  ratios  at  this  level; 
however,  if  proper  empirical  characterizations  could  be  done 
using  the  linear  model  over  many  stretch  states,  effective 
contrast  thresholds  could  be  determined  for  the  characteriza¬ 
tion  of  lesions.  In  addition,  these  results  were  also  promising 
in  that  successful  Young’s  modulus  contrast  and  localization 
was  achieved  with  a  nonpigmented  lesion.  This  indicates  that 
only  the  deflections  of  the  surrounding  image  pattern  and  not 
the  lesion  image  intensity  itself  are  responsible  for  the 
changes  in  the  elastic  modulus  values.  This  enthusiasm  must 
be  tempered  by  the  realization,  however,  that  the  in  vivo 
model  may  require  more  thought  with  respect  to  boundary 
conditions.  Undoubtedly,  the  influence  of  subcutaneous  tis¬ 
sue  connectivity  would  influence  the  results  here  if  these 
additional  constraints  were  applied.  Given  the  inherent  link 
between  the  image  formation  and  the  validity  of  the  compu¬ 
tational  model,  more  work  needs  to  be  performed  prior  to 
clinical  deployment. 
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Fig.  11.  Elasticity  image  reconstruction  of  melanoma  using  (a)  CC  and  (b)  GC  with  contrast  ratio  values  along  the  transect  for  (c)  CC  and  (d)  GC, 
respectively.  The  location  of  transect  is  designated  by  the  T  shown  in  (a)  and  (b). 


Although  these  results  are  encouraging,  not  all  reconstruc¬ 
tions  exhibited  the  same  peak  modulus  or  lesion  localization. 
One  reason  could  be  the  accuracy  to  which  boundary  condi¬ 
tions  were  determined  for  each  stretch  state.  It  is  possible 
that  the  manual  delineation  of  boundary  conditions  or  the 
observed  wrinkling  at  high  stretch  states  resulted  in  some 
boundaries  being  mapped  less  precisely  than  others.  In  some 
of  the  reconstructions,  significant  boundary  artifacts  can  be 
observed.  For  example,  in  the  second  and  fourth  column  of 


1 6  Regions  36  Regions  64  Regions 


144  Regions  256  Regions  400  Regions 


Fig.  12.  An  example  of  mutli-resolution  solution  development  using  GC  for 
melanoma  simulation. 


Fig.  7,  a  Young’s  modulus  peak  is  shown  in  the  lower  right 
hand  region  of  the  boundary.  A  second  candidate  for  recon¬ 
struction  inaccuracies  across  stretch  states  could  be  the  de¬ 
gree  of  model-data  mismatch.  It  is  interesting  to  note  the 
correlation  between  increased  stretch  and  the  marked  de¬ 
crease  in  accuracy  of  the  contrast-ratio  transect  plots.  At  the 
smaller  stretches,  3a-b  and  3a-c,  both  CC  and  GC  recon¬ 
structions  perform  better  in  both  localization  and  quantifica¬ 
tion  while  both  show  overpredictions  within  transects  for 
stretch  states,  3a-d  and  3a-e.  A  model-data  mismatch  would 
seem  a  likely  source  for  this  change  in  performance,  consid¬ 
ering  that  the  elastic  model  used  is  a  small-strain  model  and 
the  levels  of  strain  are  less  in  the  first  two  stretch  states.  One 
somewhat  qualitative  observation  that  can  also  be  made  is 
that  the  GC-based  method  appears  to  reconstruct  somewhat 
better  than  the  CC-based  method.  This  is  also  the  case  within 
the  melanoma  simulations.  Interestingly,  in  Ref.  30,  a  similar 
experience  was  found  in  that  the  GC  method  outperformed 
other  methods  with  respect  to  our  phantom  reconstructions. 
The  principal  difference  between  the  CC  and  GC  similarity 
methods  is  the  form  of  the  image  to  be  used  when  calculation 
the  correlation  coefficient.  GC  employs  the  edge  map  of  the 
image  while  CC  uses  the  raw  acquired  image.  The  increased 
performance  by  GC  may  indicate  that  areas  of  structured 
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sharp  gradient  intensities  influence  the  MIE  approach  more 
signiflcantly  than  more  gradual  intensity  changes.  However, 
this  statement  must  be  tempered  with  the  realization  of  Fig. 
10  whereby  structural  decorrelation  has  occurred  although 
arguably  at  much  larger  length  scales  as  compared  to  those  in 
traditional  USE. 

The  results  from  the  melanoma  simulations  provide  a 
more  realistic  representation  of  the  types  of  images  that  can 
be  acquired  within  the  clinic.  These  images  provide  their 
own  challenge  in  that  although  the  lesion  is  pigmented,  the 
surrounding  structured  pattern  of  the  grid  used  in  the  phan¬ 
tom  was  not  present.  In  this  case,  it  was  desirable  to  test  MIE 
without  the  presence  of  the  structured  grid.  Overall,  the  elas¬ 
ticity  images  and  transects  were  satisfying,  with  the  GC 
qualitatively  outperforming  the  CC  method.  One  interesting 
observation  is  related  to  the  apparent  suppression  of  modulus 
noise  within  the  GC  elasticity  image  as  compared  to  the  CC. 
This  is  more  than  likely  due  to  the  suppression  of  low- 
frequency  image  characteristics  associated  with  extracting 
edges  within  the  source  and  target  images. 

Figures  9  and  12  demonstrate  the  multiresolution  aspect 
to  our  approach  by  showing  the  reconstructions  at  all  six 
resolutions  used  within  the  generation  of  our  images.  One 
beneflcial  aspect  is  the  availability  of  intraresolution  elastic¬ 
ity  images  which  represent  accurate,  albeit  coarse,  assess¬ 
ments  of  image  progression.  In  addition,  these  intra¬ 
resolution  images  could  be  used  to  dynamically  alter  the 
^-means  clustering  approach  to  locally  reflne  the  reconstruc¬ 
tion  process  for  the  next  resolution  (although  not  done  in  this 
study).  This  would  alter  the  algorithm  representation  in  Fig. 
2  by  replacing  precomputed  resolution  maps  with  an  internal 
process  block  which  calculated  ^-means  regions  dynami¬ 
cally  based  on  areas  of  interest  found  during  the  reconstruc¬ 
tion  process. 

V.  CONCLUSIONS 

In  this  paper,  a  novel  multiresolution  extension  to  Modal¬ 
ity  Independent  Elastography  (MIE)  has  been  implemented 
which  simplifles  previous  work  (a  dual-grid  technique)  and 
is  shown  to  be  more  robust  than  the  single-resolution  ver¬ 
sion.  In  addition,  the  multi-resolution  architecture  imple¬ 
mented  facilitates  the  monitoring  of  reconstruction  quality  at 
intermediate  resolutions.  To  test  the  approach,  a  membrane 
experimental  setup  was  created  which  utilizes  sets  of  optical 
images  for  the  reconstruction  process.  The  use  of  optical 
images  to  generate  Young’s  modulus  reconstructions  does 
represent  a  new  modality  within  MIE  development  and  could 
potentially  be  used  within  dermoscopic  applications. 

Results  from  phantom  and  simulation  experiments  dem¬ 
onstrated  that  the  multiresolution  MIE  approach  is  viable 
within  the  context  for  both  nonpigmented  and  pigmented  le¬ 
sions,  respectively.  The  nonpigmented  phantom  experiment 
highlighted  direct  comparisons  between  images  of  Young’s 
modulus  contrast  and  their  independently  measured  counter¬ 
parts,  as  provided  by  mechanical  testing.  Overall  the  results 
indicated  good  localization  and  quantiflcation.  However,  re¬ 
sults  did  indicate  a  dependence  on  the  fldelity  of  the  recon¬ 
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struction  and  the  degree  of  applied  deformation.  In  addition 
to  the  phantom  experiment,  a  simulation  using  a  clinical  im¬ 
age  of  a  pigmented  melanoma  were  reported  and  illustrated 
excellent  localization  and  quantiflcation. 

Despite  potential  limitations  in  elasticity  image  resolution 
when  compared  to  traditional  MRE  and  USE,  MIE’s  adapt¬ 
ability  to  an  optical  image-registration  platform  at  multiple 
scales  is  an  intriguing  possibility.  Furthermore,  this  extension 
to  another  modality  demonstrates  that  MIE-based  approaches 
to  elastography  represent  a  new  class  of  algorithms  that  may 
yield  potentially  new  frameworks  for  disease  characteriza¬ 
tion. 
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